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OVERVIEW: 

The  objectives  of  this  project  are: 

1 .  Formulation  of  a  constitutive  model  for  the  analysis  of  fatigue  damage  in  metal 
matrix  composites 

2.  Development  of  a  finite  element  program  in  which  the  developed  fatigue  damage 
model  will  be  integrated 

3.  Involvement  and  training  of  graduate  students  in  this  advanced  research  in  the 
area  of  fatigue  damage  of  metal  matrix  composites 

The  appendix  contains  copies  of  scientific  articles  which  are  representative  of  the  publications 
written  during  the  course  of  this  project. 


SUMMARY  OF  ACCOMPLISHMENTS: 

The  research  accomplished  during  the  project  phase  is  described  below  in  three  major 
categories: 

I.  Theoretical  formulation  of  the  cyclic  and  fatigue  damage  model 

II.  Numerical  Implementation  of  the  model  into  a  finite  element  code 

III.  Extent  of  the  involvement  and  participation  of  graduate  students 

This  is  followed  by  a  listing  of  the  publications  written  during  the  course  of  this  project,  including 
books,  refereed  journal  articles,  conference  proceedings  and  presentations. 

I.  Theoretical  Formulation: 

The  theoretical  formulation  is  based  on  the  effective  stress  concept,  and  incorporates  the  damage- 
plasticity  theory  for  composite  materials  using  a  micromechanical  approach.  In  the  micromechanicd 
approach  the  damage  phenomenon  due  to  applied  cyclic/fatigue  loading  is  considered  at  the 
constituent  level.  For  this  purpose  a  homogenization  technique  in  the  form  of  the  Mori-Tanaka 
method  is  employed  in  order  to  allow  for  the  distribution  of  the  external  applied  loading  to  the 


individual  constituents  by  means  of  the  so-called  stress  and  strain  co'ncentration  tensors.  The 
material  behavior  is  then  modeled  at  the  constituent  level  based  on  the  distributed  loading.  Elasto- 
plastic  deformations  as  well  as  damage  initiation  and  evolution  processes  are  considered  in  the 
constituents  individually.  For  this  purpose  individual  constitutive  equations  to  model  the  elasto- 
plastic  behavior  of  and  the  damage  evolution  process  in  the  constituents  are  established  for  each 
of  the  constituents.  The  damage  evolution  process  includes  the  modelirg  of  damage  initiation  and 
damage  propagation  up  to  final  failure  for  which  a  specific  failure  criterion  is  applied.  Through  the 
use  of  the  employed  homogenization  technique  the  overall  composite  material  behavior,  damage 
state  and  damage  evolution  is  predicted  based  on  the  material  behavior  and  damage  evolution  in 
the  constituents. 

The  developed  micromechanical  cyclic/fatigue  damage  model  is  based  on  a  consistent  formulation 
in  the  framework  of  continuum  mechanics  and  on  thermodynamical  principles.  The  formulation  for 
the  damage  evolution  includes  also  the  physical  aspects  of  fatigue  damage  in  metal  matrix 
composites.  Through  the  use  of  individual  damage  variables  in  the  form  of  second  order  tensors 
for  the  constituents  it  is  possible  to  model  anisotropic  damage  development  and  anisotropic  damage 
evolution  in  the  constituents.  Furthermore,  the  micromechanical  approach  allows  to  modeldifferent 
damage  and  failure  modes  in  the  composite  material  based  on  the  failure  in  the  individual 
constituent.  Hence  composite  failure  modes  such  as  overload  failure  due  to  matrix  cracking  or  fber 
cracking,  depending  on  the  type  of  composite,  maybe  modeled.  Furthermore  tie  introduction  of  an 
additional  second  order  damage  tensor  for  the  fiber-matrix  interface  allows  for  the  modeling  of 
damage  and  failure  modes  due  to  debonding  as  well. 


II.  Numerical  Implementation: 

The  cyclic/fatigue  damage  constitutive  model  is  implemented  numerically  and  used  to  simulate 
damage  development  and  evolution  in  metal  matrix  composites.  Numerical  problems  encountered 
during  the  implementation  and  testing  phase,  such  as  numerical  instabilities,  divegence  problems, 
have  been  addressed  and  appropriate  measures  to  resolve  them  have  been  taken  and  included. 
Appropriate  necessary  material  and  model  parameters  are  established  based  on  available 
experimental  results.  The  numerical  implementation  is  used  to  simulate  various  loading  cases  for 
metal  matrix  composite  specimen,  and  the  results  are  then  compared  with  the  results  from  the 
experimental  investigations.  Comparison  of  the  results  from  the  numerical  simulations  with  those 
from  experiments  show  good  agreement  (see  sample  publications  in  the  appendix)  and 
substantiate  the  capabilities  of  the  model. 


«• 

III.  Publications,  Presentations  and  Organization  of  Symposia: 

A.  Books 

1.  Voyiadjis,  G.  Z.,  and  Allen,  D.  H.,  editors,  Studies  in  Applied  Mechanics,  Vol.  44,  Damage  and 

interfacial  Debonding  in  Composites,  275  p.,  Elsevier,  Amsterdam,  1996. 

2.  Voyiadjis,  G.  Z.,  and  Ju,  J.  W.,  editors.  Studies  in  Applied  Mechanics,  Vol.  41,  Inelasticity  and 
Micromechanics  of  Metal  Matrix  Composites,  351  p.,  Elsevier,  Amsterdam,  1994. 

3.  Voyiadjis,  G.  Z.,  Bank,  L.  C.,  and  Jacobs,  L.  J.,  editors.  Studies  in  Applied  Mechanics,  Vol.  35, 
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review). 
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12.  Voyiadjis,  G.  Z.,  and  Thiagarajan,  G.,  "A  Cyclic  Plasticity  Model  for  Metal  Matrix  Composites 
Using  an  Anisotropic  Yield  Surface,”  Studies  in  Applied  Mechanics,  Vol.  41,  Inelasticity  and 
Micromechanics  of  Metal  Matrix  Composites,  edited  by  G.  Z.  Voyiadjis  and  J.  W.  Ju,  Elsevier 
Publishing  Company,  1994,  pp.51-81. 

13.  Voyiadjis,  G.  Z.,  and  Park,  T.,  "Interfacial  and  Local  Damage  Analysis  in  Metal  Matrix 
Composites,"  Damage  Mechanics  in  Composites,  1994,  edited  by  D.H.  Allen,  and  J.W.  Ju, 

ASME,  AMD-Vol.  185,  1994,  pp.  87-108. 

14.  Voyiadjis,  G.  Z.,  and  Kattan,  P.  I.,  "Coupling  of  Damage  and  Viscoplasticity  for  Large  Deformation 
of  Metals,"  Large  Plastic  Deformations;  Fundamental  Aspects  and  Applications  to  Metal  Forming, 
edited  by  C.  Teodosiu,  J.  L.  Raphanel,  and  F.  Sidoroff,  A.  A.  Balkema  Publishing  Company,  1993, 
pp.  345-352. 

15.  Voyiadjis,  G.  Z.,  Kattan,  P.  I.,  and  Venson,  A.  R.,  "Evaluation  of  a  Damage  Tensor  for  Metal 
Matrix  Composites,"  Direction  des  Etudes  et  Recherches  d’  Electricite  de  France;  MECAMAT  93, 
International  Seminar  on  Micromechanics  of  Materials,  Vol.  84, 1993,  pp.  406-417. 

16.  Voyiadjis,  G.  Z.,  and  Kattan,  P.  I.,  "Micromechanical  Characterization  of  Damage-Plasticity  in 
Metal  Matrix  Composites,"  Studies  in  Applied  Mechanics,  Vol.  34,  Damage  in  Composite 
Materials,  edited  by  G.  Z.  Voyiadjis,  Elsevier  Publishing  Company,  1993,  pp.  67-102. 


D.  Technical  Reports 


1.  Voyiadjis,  G.  2.,  Kattan,  P.  I.,  Venson,  A.  R.,  Park,  T.,  and  Thiagarajan,  G.,  "Constitutive  Modeling 
of  Metal  Matrix  Composites  Under  Cyclic  Loads  -  Second  Annual  Report  (1993-94),  F49620-93-1- 
•  0097DEF  and  F49620-92-J-0463,"  Air  Force  Office  of  Scientific  Research,  Washington,  D.C., 

September  1994, 155  pages. 
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• 

2.  Voyiadjis,  G.  Z.,  Kattan,  P.  I.,  and  Venson,  A.  R.,  "Constitutive  Modeling  and  Finite  Element 

Analysis  of  Fatigue  Damage  in  Metal  Matrix  Composites  -  First  Annual  Report  (1992-93),  F49620- 
93-1-0097DEF  and  F49620-92-J-0463,"  Air  Force  Office  of  Scientific  Research,  Washington, 

D.C.,  October  1993, 130  pages. 

•  3.  Voyiadjis,  G.  Z.,  Kattan,  P.  I.,  Venson,  A.  R.,  and  Park,  T.,  "Anisotropic  Damage  Mechanics 

Modeling  in  Metal  Matrix  Composites  -  Final  Report,  AFOSR-90-0227  DEF,"  Air  Force  Office  of 
Scientific  Research,  Washington,  D.C.,  May  1993, 141  pages. 

4.  Voyiadjis,  G.  Z.,  Kattan,  P.  I.,  and  Venson,  A.,  "Anisotropic  Damage  Mechanics  Modeling  in  Metal 
Matrix  Composites  -  Second  Annual  Report  (1991-92),  AFOSR-90-0227  DEF,"  Air  Force  Office  of 
^  Scientific  Research,  Washington,  D.C.,  May  1992, 104  pages. 


E.  Papers  presented  in  conferences  and  invited  lectures 

1.  Voyiadjis,  G.  Z.,  and  Echle,  R.,  "Fatigue  Damage  in  Metal  Matrix  Composites."  Invited  lecture 
presented  in  the  "Symposium  on  Failure  Mechanisms  and  Mechanisms-Based  Modeling  in  High 
Temperature  Composites,"  of  the  1996  ASME,  International  Mechanical  Engineering  Congress 
and  Exposition,  Atlanta,  Georgia,  November  1996. 

2.  Voyiadjis,  G.  Z.,,  and  Echle,  R.,  "A  Micro-Mechanical  Fatigue  Damage  Model  for  Uni-Directional 
Metal  Matrix  Composites."lnvited  lecture  presented  in  the  "High  Temperature  Fatigue  in  MMCs 
Symposium,"  of  the  33rd  Annual  Technical  Meeting  of  the  Society  of  Engineering  Science, 
Arizona  State  University,  Tempe,  Arizona,  October  1996. 

3.  Voyiadjis,  G.  Z.,  and  Park,  T.,  "Quantification  of  Different  Types  of  Damage  in  Constituents  of 
MMC  Laminates."  Invited  lecture  presented  at  the  "Prager  Symposium"  of  the  33rd  Annual 
Technical  Meeting  of  the  Society  of  Engineering  Science,  Arizona  State  University,  Tempe, 
Arizona,  October  1996. 

4.  Voyiadjis,  G.  Z.,  and  Park,  T.,  "The  Kinematics  of  Damage  in  Metal  Matrix  Composites."  Invited 
lecture  presented  at  the  First  Euroconference  and  U.S.  Workshop  on  Material  Instabilities  in 
Deformation  and  Fracture,  Porto  Carras,  Chalkidiki,  Greece,  September  1996. 

5.  Echle,  R.,  and  Voyiadjis,  G.  Z.,  "A  Micro-Mechanical  Fatigue  Damage  Approach  for  Metal  Matrix 
Composites."  Presented  at  the  "Third  International  Conference  on  Composites  Engineering, 
ICCE/3,  New  Orleans,  Louisiana,  July  1996. 

6.  Voyiadjis,  G.  Z.,  and  Deliktas,  6.,  "Local  Damage  using  the  Generalized  Cell  Method  in  MMCs." 
Invited  lecture  presented  at  the  "Third  International  Conference  on  Composites  Engineering," 
ICCE/3,  New  Orleans,  Louisiana,  July  1996. 

7.  Park,  T.,  and  Voyiadjis,  G.  Z.,  "Interfacial  Damage  in  Metal  Matrix  Composites."  Presented  at  the 
'Third  International  Conference  on  Composites  Engineering,"  ICCE/3,  New  Orleans,  Louisiana, 
July  1996. 

8.  Thiagarajan,  G.,  and  Voyiadjis,  G.  Z.,  "In-Situ  Modeling  of  Titanium  Based  Matrix  in  Metal-Matrix 
Composites  Subjected  to  Cyclic  Loading."  Presented  at  the  "Third  International  Conference  on 
Composites  Engineering,"  ICCE/3,  New  Orleans,  Louisiana,  July  1996. 
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9. 


10. 


11. 


12. 


13. 


14. 


15. 


16. 


•  17. 


18. 


19. 


Voyiadjis,  G.  Z.,  and  Deliktas,  B.,  "Damage  in  Metal  Matrix  Composites  Using  the  Generalized 
Cells  Method."  Presented  at  the  1996  ASME  Mechanics  and  Materials  Conference,  Baltimore, 
Maryland,  June  1996. 

Voyiadjis,  G.  Z.,  and  Echle,  R..  "Fatigue  Model  for  Damage  in  Metal  Matrix  Composites."  Invited 
lecture  presented  in  the  Symposium  on  "Application  of  Continuum  Damage  Mechanics  to  Fatigue 
and  Fracture,"  at  the  ASTM  Conference,  Orlando,  Florida,  May  1996. 

Voyiadjis,  G.  Z.,  and  Park,  T.,  "Effect  of  Damage  Due  to  Delamination  on  the  Flexural  Stiffness  of 
Laminated  Composite  Plates."  Invited  lecture  presented  at  the  Symposium  on  "Damage  and 
Interfacial  Debonding  in  Composites,"  at  the  32nd  Annual  Technical  Meeting  of  the  Society  of 
Engineering  Science,  Ne\w  Orleans,  Louisiana,  October  1995. 

Thiagarajan,  G.,  and  Voyiadjis,  G.  Z.,  "Cyclic  Damage  with  Cyclic  Plasticity  in  Metal  Matrix 
Composites."  Lecture  presented  at  the  32nd  Annual  Technical  Meeting  of  the  Society  of 
Engineering  Science,  New  Orleans,  Louisiana,  October  1995. 

Voyiadjis,  G.  Z.,  and  Guelzim,  Z.,  "Incremental  Damage  Theory  for  Metal  Matrix  Composites." 
Invited  lecture  presented  at  the  "Eringen  Medal  Symposium:  S.  Atiuri,"  at  the  32nd  Annual 
Technical  Meeting  of  the  Society  of  Engineering  Science,  New  Orleans,  Louisiana,  October  1995. 

Park,  T.,  and  Voyiadjis,  G.  Z.,  "Analysis  of  Interfacial  and  Local  Damage  in  Metal  Matrix 
Composites  using  the  FE  Method."  Lecture  presented  at  the  32nd  Annual  Technical  Meeting  of 
the  Society  of  Engineering  Science,  New  Orleans,  Louisiana,  October  1995. 

Thiagarajan,  G.,  and  Voyiadjis,  G.  Z.,  "Cyclic  Plasticity  with  Damage  in  Metal  Matrix  Composites." 
Lecture  presented  at  the  Symposium  on  "Micromechanics  of  Fibers  or  Composites,"  at  the 
Second  International  Conference  on  Composites  Engineering,  New  Orleans,  Louisiana,  August 
1995. 

Voyiadjis,  G.  Z.,  and  Rainer,  E.,  "Fatigue  Damage  in  Metal  Matrix  Composites."  Invited  lecture 
presented  at  the  Symposium  on  "Fatigue  of  Composites,"  at  the  Second  International  Conference 
on  Composites  Engineering,  New  Orleans,  Louisiana,  August  1995. 

Park,  T.,  and  Voyiadjis,  G.  Z.,  "Stress  and  Strain  Concentration  Tensors  for  Interfacial  Damage." 
Lecture  presented  at  the  Symposium  on  "Micromechanics  of  Fibers  or  Composites,"  at  the 
Second  International  Conference  on  Composites  Engineering,  New  Orleans,  Louisiana,  August 
1995. 

Venson,  A.  R.,  and  Voyiadjis,  G.  Z.,  "Experimental  Determination  of  Damage  in  Titanium 
Aluminide  with  Silicon  Carbide  Fibers."  Lecture  presented  at  the  Symposium  on  "Metal  Matrix 
Composites,"  at  the  Second  International  Conference  on  Composites  Engineering,  New  Orleans, 
Louisiana,  August  1995. 

Voyiadjis,  G.  Z.,  and  Thiagarajan,  G.,  "A  Cyclic  Plasticity  Model  for  Metal  Matrix  Composites." 
Invited  lecture  presented  at  the  Symposium  on  "Recent  Advances  in  Constitutive  Laws  for 
Engineering  Materials,"  at  the  International  Conference  on  Computational  Engineering  Science 
{ICES-95),  Mauna  Lani,  Hawaii,  August  1995. 
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Voyiadjis,  G.  Z.,  and  Park,  T.,  "Finite  Element  Analysis  of  Damage  and  Elasto-Plastic  Behavior  of 
Metal  Matrix  Composites."  Invited  lecture  presented  at  the  Symposium  on  "Inelastic  Deformation, 
Damage  and  Life  Analysis,"  at  the  International  Conference  on  Computational  Engineering 
Science  (ICES-95),  Mauna  Lani,  Hawaii,  August  1995. 

Voyiadjis,  G.  Z.,  and  Thiagarajan,  G.,  "A  Damage  Cyclic  Plasticity  Model  for  Metal  Matrix 
Composites."  Invited  lecture  presented  at  the  Symposium  on  "Micromechanics  and  Constitutive 
Modeling  of  Composite  Materials,"  at  the  ASME  joint  Applied  Mechanics  and  Materials  Summer 
Conference  AMD-MD  '95,  Los  Angeles,  California,  June  1995. 

Voyiadjis,  G.  Z.,  and  Park,  T.,  "Finite  Element  Analysis  of  Interfacial  and  Local  Damage  in  Metal 
Matrix  Composites."  Invited  lecture  presented  at  the  Symposium  on  "Recent  Advances  in 
Damage  Mechanics,"  at  the  ASME  joint  Applied  Mechanics  and  Materials  Summer  Conference 
AMD-MD  '95,  Los  Angeles,  California,  June  1995. 

Voyiadjis,  G.  Z.,  Park,  T.,  and  Venson,  A.  R.,  "Finite  Element  Implementation  of  the  Overall 
Approach  to  Damage  in  Metal  Matrix  Composites."  Invited  lecture  presented  at  the  Symposium  on 
"Damage  Mechanics  in  Engineering  Materials,"  at  the  10th  ASCE  Engineering  Mechanics 
Conference,  Boulder,  Colorado,  May  1995. 

Voyiadjis,  G.  Z.,  "Damage  in  Materials."  Invited  lecture  presented  at  Louisiana  State  University, 
Saturday  Science  at  LSU  Seminar  Series,  Baton  Rouge,  Louisiana,  December  1994. 

Voyiadjis,  G.  Z.,  and  Park,  T.,  "Interfacial  and  Local  Damage  Analysis  of  Metal  Matrix 
Composites."  Invited  lecture  presented  at  the  Symposium  on  "Damage  Mechanics  in 
Composites,"  at  the  ME  '94,  The  International  Mechanical  Engineering  Congress  and  Exposition, 
Chicago,  Illinois,  November  1994. 

Voyiadjis,  G.  Z.,  "Constitutive  Modeling  of  Plasticity  with  Damage  in  Metal  Matrix  Composites." 
Invited  lecture  presented  at  the  National  Science  Foundation  Workshop  on  Mechanics  and 
Processing  of  Advanced  Engineering  Materials  -  II,  Atlanta,  Georgia,  October,  1994. 

Voyiadjis,  G.  Z.,  Venson,  A.R.,  and  Kattan,  P.l.  "Damage  Approaches  to  Unidirectional 
Elastoplastic  Laminas  Under  Uniaxial  Tension."  Invited  lecture  presented  at  the  31st  Annual 
Technical  Meeting  of  the  Society  of  Engineering  Science,  College  Station,  Texas,  October  1994. 

Voyiadjis,  G.  Z.,  and  Kattan,  P.  I.,  "Damage  and  Plastic  Deformation  in  a  Uniaxially  Loaded  Thin 
Composite  Ply."  Invited  lecture  presented  at  the  First  International  Conference  on  Composites 
Engineering,  New  Orleans,  Louisiana,  August  1994. 

Voyiadjis,  G.  Z.,  and  Kattan,  P.  I.,  "Micromechanical  Modeling  of  Damage  and  Plasticity  in 
Continuously  Reinforced  Metal  Matrix  Composites."  Invited  lecture  presented  at  the  12th  U.S. 
National  Congress  of  Applied  Mechanics,  Seattle,  Washington,  June  1994. 

Voyiadjis,  G.  Z.,  "Damage  Accumulation  in  Metal  Matrix  Composites."  Invited  lecture  presented  at 
the  National  Science  Foundation  Workshop  on  "Processing  and  Constitutive  Modeling  of 
Advanced  Engineered  Materials  - 1,"  Washington  State  University,  Pullman,  Washington,  October 
1993. 


# 


31 .  Voyiadjis,  G,  Z.,  "Damage  in  Metal  Matrix  Composites."  Invited  lecture  presented  in  the 

Department  of  Mechanical  Engineering,  University  of  Delaware,  Newark,  Delaware,  October 
1993. 


•  32.  Voyiadjis,  G.  Z.,  Venson,  A.  R.,  and  Kattan,  P.  I.,  "Damage  Mechanism  in  Metal  Matrix  Composite 

Plates."  Invited  lecture  presented  at  the  Army  Research  Office  Workshop  on  Dynamic  Response 
of  Composite  Structures,  New  Orleans,  Louisiana,  August  1993. 

33.  Voyiadjis,  G.  Z.,  Kattan,  P.  I.,  and  Venson,  A.  R.,  "Evaluation  of  a  Damage  Tensor  for  Metal 

Matrix  Composites."  Presented  at  the  MECAMAT  '93,  International  Seminar  on  Micromechanics 
of  Materials,  Fontainebleau,  France,  July  1993. 


34.  Voyiadjis,  G.  Z.,  and  Sivakumar,  S.  M.,  "Cyclic  Plasticity  and  Ratchetting."  Presented  at  the 
MEErN'93,  First  SES-ASME-ASCE  Joint  Meeting  on  Mechanics,  Bieniek  Symposium  on 
Mechanics  of  Materials  and  Structures,  Charlottesville,  VA,  June  1993. 


35.  Voyiadjis,  G.  Z.,  and  Kattan,  P.  I.,  "Damage  and  Inelastic  Deformation  of  SiC/Ti-Aluminide 
®  Composites."  Invited  lecture  presented  at  the  MEET'N'93,  First  SES-ASME-ASCE  Joint  Meeting 

on  Mechanics,  Symposium  on  Inelastic  Micromechanics  in  SiC/Ti  Composites,  Charlottesville, 
VA,  June  1993. 


F. 


1. 


2. 


3. 


4. 


5. 


6. 


7. 


Symposia  organized  by  Dr.  G.  Z.  Voyiadjis: 

Co-Organizer  and  Co-Chairman  of  three  sessions  of  the  Symposium  on  "Failure  Predictions  in 
Dynamic  Environments,"  to  be  held  in  the  1997  ASME,  IMECE,  Dallas,  Texas,  November  1997. 

Co-Organizer  and  Co-Chairman  of  eleven  sessions  of  the  Symposium  on  "Damage  Mechanics  in 
Engineering  Materials,"  to  be  held  in  the  Joint  ASME/ASCE/SES  Mechanics  Conference, 
Northwestern  University,  Evanston,  Illinois,  June/July  1997. 

Organizer  and  Chairman  of  the  Symposium  on  "High  Temperature  Fatigue  in  Metal  Matrix 
Composites,"  held  in  the  Society  of  Engineering  Science  33rd  Annual  Technical  Meeting,  Arizona 
State  University,  Tempe,  Arizona,  October  1996. 

Co-Organizer  and  Co-Chairman  of  four  sessions  of  the  Symposium  on  "Damage  and  Interfacial 
Debonding  in  Composites,"  held  in  the  32nd  Society  of  Engineering  Science  Meeting,  New 
Orleans,  Louisiana,  October/November  1995. 

On  the  Organizing  Committee  of  the  32nd  Annual  Technical  Meeting  of  the  Society  of 
Engineering  Science,  New  Orleans,  Louisiana,  October  29  -  November  1, 1995. 

Co-Organizer  and  Co-Chairman  of  four  sessions  of  the  Symposium  on  "Micromechanics  of  Fibers 
or  Composites,"  held  in  the  2nd  International  Conference  on  Composites  Engineering,  New 
Orleans,  Louisiana,  August  1995. 

Co-Organizer  and  Co-Chairman  of  four  sessions  of  the  Symposium  on  "Fatigue  of  Composites," 
to  be  held  in  the  2nd  International  Conference  on  Composites  Engineering,  New  Orleans, 
Louisiana,  August  1995. 
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8.  Co-Organizer  and  Co-Chairman  of  five  sessions  of  the  Symposium  on  "Inelasticity  and 
Micromechanics  of  Metal  Matrix  Composites,"  held  in  the  Twelfth  U.S.  National  Congress  of 
Applied  Mechanics,  Seattle,  Washington,  June/July  1994. 

9.  Panel  Member  of  the  session  on  "Damage  Accumulation  and  Life  Prediction  in  Metal  Matrix 
Composites,"  of  the  National  Science  Foundation  Workshop  on  "Processing  and  Constitutive 
Modeling  of  Advanced  Engineered  Materials,"  Washington  State  University,  Pullman, 
Washington,  October  1993. 

10.  Co-Organizer  and  Co-Chairman  of  four  sessions  of  the  Maciei  P.  Bieniek  Symposium  on 
"Mechanics  of  Materials  and  Structures,"  held  in  the  MEET'N'93,  First  Joint  ASCE-ASME-SES 
Mechanics  Conference,  University  of  Virginia,  Charlottesville,  Virginia,  June  1993. 


Ph.D.  Students  Completed: 

1 .  Rainer  Echle:  A  Micro  -  Mechanical  Fatigue  Damage  Model  for  Uni  -  Directional  Metal  Matrix 
Composites;  May  1997,  LSU. 

2.  Ganesh  Thiagarajan:  A  Cyclic  Plasticity/Damage  Model  for  Metal  Matrix  Composites;  May 
1996,  LSU.  Currently:  Instructor,  Department  of  Civil  and  Environmental  Engineering, 
Louisiana  State  University,  Baton  Rouge,  Louisiana. 

3.  A.  R.  Venson:  Experimental  Macro  and  Microstructural  Characterization  of  Damage  for 
Metal  Matrix  Composites;  December  1994,  LSU.  Currently:  Assistant  Professor,  Department 
of  Civil  Engineering,  University  of  Southwestern  Louisiana,  Lafayette,  Louisiana. 


4. 


Taehyo  Park:  Finite  Element  Analysis  of  Damage  and  Elastic-Plastic  Behavior  of  Metal 
Matrix  Composites:  December  1994,  LSU.  Currently:  Research  Associate,  Department  of 
Civil  and  Environmental  Engineering,  Louisiana  State  University,  Baton  Rouge,  Louisiana 


M.Sc.  Students  Completed: 


Babur  Deliktas:  Damage  in  Metal  Matrix  Composite  Using  the  Generalized  Cells  Model;  May 
1996,  LSU. 
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The  following  copies  of  scientific  papers  are  representative  of  the  publications  written  during  the 
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A  CYCLIC  ANISOTROPIC-PLASTICITY  MODEL  FOR  METAL 

MATRIX  COMPOSITES 

George  Z.  Voyiadjis  and  Ganesh  Thiagarajan 
Department  of  Civil  Engineering,  Louisiana  State  University,  Baton  Rouge,  LA-70803,  U.S.A. 


{Received  in  final  revised  form  5  January  1995) 


Abstract— Based  on  a  six  parameter  general  anisotropic  yield  surface  proposed  earlier  by  Voyiadjis 
and  Thiagarajan  (An  Anisotropic  Yield  Surface  Model  for  Directionally  Reinforc^  Metal  Matrix 
Composites,  Int.  J.  Plasticity  [1995]),  a  cyclic  plasticity  model  to  model  the  behavior  of  direction- 
ally  reinforced  metal  matrix  composite,  has  been  proposed  here.  Apart  from  being  able  to  model 
®  different  initial  yielding  behavior  along  different  stress  directions,  a  number  of  features  have  been 

incorporated  into  the  plasticity  model.  They  include  the  usage  of  a  proposed  non-associative  flow 
rule,  kinematic  hardening  rule  of  Phillips  type,  a  modified  form  of  the  bounding  surface  model  for 
modelling  the  cyclic  behavior,  and  the  usage  of  a  proposed  form  for  evaluating  the  plastic  modtflus 
for  anisotropic  materials.  Previous  experimental  data  have  been  used  for  the  evaluation  of  the  yield 
surface  parameters  as  well  as  those  for  the  determination  of  the  plastic  modulus.  The  stress— strain 
_  resulu  generated  from  the  model  have  then  been  compared  with  those  from  the  experiments.  The 

•  behavior  of  the  model  under  certain  simulated  cyclic  loading  situations  has  also  been  presented. 

1.  INTRODUCTION 

This  paper  focuses  upon  the  treatment  of  a  metal  matrix  composite  (MMC)  as  a  con- 

•  tinuum  and  is  an  extension  of  the  work  presented  earlier  by  the  authors  (Vojiadjis  & 

Thiagarajan  [1995])  wherein  an  anisotropic  yield  surface  has  been  proposed  wth 
application  to  continuous  directionally  reinforced  metal  matrix  composites.  The  yield 
surface  proposed  earlier  has  been  correlated  to  the  experimental  observations  of 
Dvorak  et  al.  [1988]  and  Nigam  et  al.  [1993].  A  model  is  proposed  herein  to  account  for 

•  the  material  behavior  of  MMC  which  is  transversely  isotropic  and  subjected  to  cyclic, 

proportional  and  non-proportional  loadings.  The  cyclic  plasticity  model  for  the  aniso¬ 
tropic  material  is  based  on  a  modification  of  the  bounding  surface  model  proposed  by 
Dafalias  and  Popov  [1976]  for  the  case  of  isotropic  materials. 

It  is  also  observed  that  the  plastic  strains  that  develop  in  an  MMC  are  non-associative 

•  in  nature.  To  account  for  this,  a  non-associative  flow  rule  is  proposed  here,  based  on 

definitions  of  the  yield  and  a  complimentary  yield  function,  termed  as  “constrained 
yield  function.”  Also  a  suitable  kinematic  hardening  rule  is  adopted  here. 

Based  on  the  above  observations  the  elastoplastic  stiffness  matrix  is  derived  for  the 
cyclic  plasticity  model.  This  model  has  been  implemented  in  a  computer  program  to 
^  generate  the  stress-strain  curves  under  different  loading  conditions.  These  curves  have 

been  generated  along  loading  paths  followed  by  Nigam  et  al.  [1993].  Plastic  strains  are 
obtained  and  compared  to  those  documented  by  Dvorak  et  al.  [1988]  and  Nigam  et  al. 
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[1993].  Other  simulated  loading  conditions  are  applied  to  the  proposed  model  and 
results  are  presented  here,  to  demonstrate  the  versatility  of  the  model  under  various 
loading  conditions. 

n.  PROPOSED  YIELD  SURFACE 

Voyiadjis  and  Thiagarajan  [1995]  proposed  an  anisotropic  yield  surface  for  con¬ 
tinuous  fiber  reinforced  metal  matrix  composites,  using  a  fourth  order  anisotropic  yield 
tensor  M.  Two  coordinate  systems  are  used,  namely  the  local  coordinate  system  and 
the  global  coordinate  system  (Fig.  1).  In  the  local  coordinate  system,  the  coordinate 
axes  coincide  with  the  three  principal  axes  of  material  anisotropy.  The  orientation  of 
the  fiber  is  defined  in  the  global  coordinate  system.  The  loading  is  also  defined  with 
respect  to  the  global  coordinate  system. 

The  proposed  yield  function  in  the  local  coordinate  axes  is  of  the  form  shown  below. 

Mijki&ijffki  “1=0  (1) 

where  Cij  is  the  overall  state  of  stress  in  this  system.  M  is  the  fourth  order  anisotropic 
yield  tensor  expressed  as  a  function  of  two  second  order  tensors  ay  and  by  as  follows: 

M  =  M(a,b)  (2) 

The  expression  for  M  is  given  as, 

Mijki  =  AiflijOki)  +  B{aikaji)  +  C{ayajk)  -f-  D^bybu)  (3) 

where  A,  C  and  D  are  constants  and  ay  and  by  are  functions  of  the  six  strength 
parameters  kj  (/=!,  ...,  6).  Three  of  these  parameters  are  directly  related  to  the  axial 
strengths  and  the  other  three  are  shear  strength  parameters.  They  are  used  to  define 
yielding  for  an  anisotropic  material.  These  parameters  are  measured  and  determined 
along  the  principal  axes  of  anisotropy,  ay  and  by  are  given  as  follows: 

>1  0  O' 

Oij  =  0  ^2  0  (4) 

0  0  ifeaj 

'0  kA  ks 

by  —  ki  0  .  (5) 

^5^6  0 

Substituting  (3)  into  (1)  we  can  express  the  yield  equation  in  the  local  coordinate 
axes  in  component  form  as  foUows: 

{A-h  B  +  C)(A:j^i  -I-  +  ^^3^3) 

+{2A){kik2&na22  +  ^1^3^11^33  +  ^2^3^22^33) 

C)kik2  +  ADki)d]2  +  (2(B+  C)kik3  +  4i)A|)^3 

+i2{B  +  C)k2h  +  =  0. 


(6) 


The  constants  A,  B,  C  and  D  are  not  material  parameters.  The  possible  choices  of 
combinations  of  these  values  are  outlined  by  Voyiadjis  and  Thiagarajan  n99^.  The 
values  of  these  constants  chosen  here  for  this  implementation  zrt  B-C-^  and 

Z)=4.  These  values  reduce  the  above  equation  to  the  following  form. 

^  =  5  (*1^1  +  ^2^22  +  ^^3) 

—  ^(^1^2^11^22  +  ^2^3^22^33  +^1^3^11^33) 

+ 1  {k\k2  +  +  I  (^1*3  +  *^)^3  +  5  (^2^3  +  -  l-O- 

One  can  also  express  the  yield  equation  in  the  global  axes  of  reference  as  follows: 

OijMijkl<^kl  “1=0.  (8) 

The  stresses  are  transformed  from  the  global  to  the  local  axes  of  reference  as  follows 

Oij  =  djpffpqdqj  (^) 

where  dy  are  the  coefficients  of  the  orthogonal  transformation  matrix.  Assuming  that 
the  fibers  are  aligned  along  the  j:-axis  (1 -direction)  one  can  write, 

dij  =  {vi,m^m) 

where  t;,,  (i  =  1,2,3)  are  the  direction  cosines  of  the  fiber  in  the  global  coordinate  system. 
A  lamina  of  any  arbitrary  orientation  is  derived  by  rotating  the  principal  axes  of  am- 
sotropy  about  the  z-axis(3-direction)  of  the  global  axes  of  reference.  Hence  one  obtains 
dy’^iO,  0,1).  The  terms  of  ^^2;  can  then  be  derived  from  the  condition 

dpid,i  =  6p,.  (11) 

Substituting  for  ay  in  the  yield  (1)  one  obtains 

pqdipdjq^d ij}^ldi^fndln^fnn  “  1  ■“  0.  (1^) 

From  the  above  equation  Miju  can  be  expressed  as, 

Mijkl  —  Mpqrs  dipdjqdkrdis-  (13) 

Further  details  about  the  computation  of  the  parameters,  their  reduction  to  certain  well- 
known  yield  criteria  and  their  comparison  with  the  yield  surfaces  obtained  from  experi¬ 
mental  evidence  have  been  demonstrated  in  the  paper  by  Voyiadjis  and  Thiagarajan  [1995]. 

m.  THEORETICAL  DEVELOPMENT  OF  THE  CONSTITUTIVE  MODEL 

The  description  of  the  elasto-plastic  behavior  of  the  metal  matrix  composite,  when 
treated  as  a  continuum  is  a  complex  task.  The  fact  that  MMC  is  transversely  isotropic 
and  the  presence  of  continuous  fibers  with  their  respective  constraints,  necessitates  the 
usage  of  anisotropic  hardening  and  non-associative  flow  rules.  The  formulation  pro¬ 
posed  here  is  also  intended  to  describe  the  behavior  of  MMC  under  cyclic  loading 
conditions.  The  bounding  surface  model  is  adopted  here  to  simulate  the  mechanical 
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behavior  of  the  material.  In  this  paper,  the  thrust  is  thus  on  transversely  isotropic 
materials. 


III.  1 .  Elastic  behavior 

The  elastic  behavior  of  the  composite  material,  treated  as  a  homogeneous  continuum 
with  transversely  isotropic  properties  has  been  defined  Walpole  [1969]  and  is  used  here. 
The  linear  constitutive  relation  is  expressed  as 

where  C  is  the  fourth  order  elastic  stiffness  tensor  relating  the  symmetric  second  order 
tensors  c  and  c  of  stress  and  strain,  respectively.  For  a  transversely  isotropic  material 
the  fourth  order  elastic  stiffness  tensor  is  given  as  follows: 


Cijki  —  Ktijtki  +  +  2pE*ji^ 

where, 

tij  =  ntij  +  2l//|y 

A;  =  ViW 

ntij  —  ^ij  ~ 

^]jki  -  \  [fnikniji  +  mjknta  -  mym*;] 

^*jki  =  \  [^ikhi  +  ^iihk  +  +  nijkma] 


(15)  ^ 

(16) 

(17) 

(18)  • 

(19) 

(20) 


and  ‘K’  is  the  plane>strain  bulk  modulus,  w,  is  the  transverse  shear  modulus,  p  is  the 
axial  shear  modulus  and  E  and  u  are  the  Young’s  Modulus  and  Poisson  ratio,  respec¬ 
tively,  when  the  material  is  loaded  in  the  fiber  direction.  For  a  transversely  isotropic 
material  the  plane-strain  bulk  modulus  can  be  defined  in  terms  of  the  other  four  elastic 
constants. 


III.2.-  Non-associative  flow  rule 

It  has  l^n  well  established  by  now  that  the  determination  of  plastic  strains,  for  any 
anisotropic  material  in  general,  and  a  MMC  in  particular  must  adopt  a  non-associative 
flow  rule.  Dvorak  et  al.  [1988]  and  Nigam  et  al.  [1993]  have  experimentally  demonstrated 
this  and  have  established  that  there  is  a  tendency  for  the  direction  of  plastic  strains  to 
be  more  inclined  towards  the  shear  direction  in  a  combined  transverse  tension-shear 
loading  situation.  It  has  also  been  shown  that  there  exists,  to  a  large  extent,  plastic 
inextensibility  along  the  fiber  direction. 

To  account  for  the  above  factors,  a  plastic  potential  function  is  proposed  here,  the 
form  of  which  is  based  on  the  proposed  yield  function.  The  nature  of  this  potential 
function  is  explained  below. 


III.3.  Proposed  potential  function 

In  order  to  determine  the  plastic  strain  increments  (c|))  use  is  made  of  a  non- 
associative  flow  rule  as  outlined  earlier  such  that. 
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dffij 


(21) 


where  G  is  the  plastic  potential  function.  The  potential  function  is  defined  later  as  a 
function  of  the  yield  function  and  the  constrained  yield  function.  The  consttained  yield 
function  is  defined  such  that  it  satisfies  the  condition  of  plastic  inextensibility  along  the 
direction  defined  by  r\.  This  is  accomplished  by  defining  a  function  g  which  is  of  the 
functional  form  of  the  yield  function/.  The  function  g  is  defined  usmg  the  fourth  order 
anisotropic  yield  trasor  M  and  a  constrained  stress  term  such  that, 

g  =  rijMijkiru  -  1  •  (^^) 


The  constraint  that  is  introduced  in  the  stress  term  is  that  the  plastic  stram  increment  is 
independent  of  the  component  of  stress  along  a  specified  direction  (defined  here  by »?/). 
Following  the  procedure  outlined  by  Spencer  [1972]  the  constraint  is  incorporated  into 
the  stress  term  as  follows: 

r,7  =  (Tij  -  Trjirjj  (23) 


where  TTjiTu  is  the  reaction  to  an  inextensibility  constraint  along  the  direction  r\.  Taking 
the  inner  product  on  both  sides  of  (23)  with  rjiTij  one  obtains, 

ri/mVj  =  “  TriJifniiy.  (24) 

One  can  impose  the  constraint  is  in  the  stress  term  ry  as  follows 

nfitrij  =  0.  (25) 


Substituting  this  constraint  in  (24)  one  obtains, 

=  Triiiyrnry  (26) 

and  using  the  condition  rjiiy  =  1  it  can  be  shown  that 

r=(T//T?,TJ.  (27) 


Hence  one  can  expand  (23)  as  follows: 

fij  =  (Tij-  {arznrVs)riirij-  (2^) 

Using  the  yield  function  and  the  function  g  defined  above,  the  potential  function  G 
c-?n  now  be  defined  as  follows. 


•  G  =  u;/+(1.0-w)g,  0<a;<1.0. 


(29) 


In  order  to  illustrate  the  capability  of  this  potential  function  we  consider  the  yield 
function  in  the  an  -  an  space.  Assuming  that  rj  -  (1,0,0)  then /and  g  functions  can  ^ 
represented  as  shown  in  Fig.  2.  Since  t)  *=(1,0,0),  this  implies  that  the  fiber  is  along  the 
1 -direction  and  the  stress  along  this  direction  an  does  not  influence  yielding.  Hence  m 
Fig.  2  the  yield  function /represents  yielding  both  axes  whereas  the  constrained  yield 
function  g  is  parallel  to  the  an-axis,  physically  representing  the  requirement  that  an 
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does  not  influence  yielding.  The  unit  normals  to  these  functions,  which  are  also  shown 
in  Fig.  2,  can  be  expressed  in  tensor  form  as  follows: 


x-y:  general  oxes  of  reference 
1-2:  principol  moteriol  axes 
- direction  of  fiber 

Fig.  1.  Local  and  global  axes  of  reference  for  a  single  lamina. 


Yield  Surface  T 

—  Constrained  Yield  Surface 'g* 
with  Ti  *(1,0,0) 

Fig.  2.  Illustration  of  non-associative  flow  rule  in  (tji  -  space. 
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The  non-associativity  of  the  flow  rule  in  built  into  the  definition  of  the  potential 
function  through  the  factor  u.  Based  on  the  flow  rule,  one  can  define  a  second  order 
tensor  representing  the  direction  of  plastic  strains  as  follows. 

mj  =  (1 .0  - u})n]j  .  (32) 

The  experimental  evidence  that  the  plastic  strains  are  predominantly  along  the  shear 
direction  can  be  incorporated  into  the  flow  rule  by  using  a  value  of  u;  in  between  0.0  and 

#  1.0.  A  value  of  a;  =  1.0  gives  a  purely  associative  flow  rule.  Using  a  value  of  w— 0.0  would 

result  in  the  usage  of  only  the  function  g  in  (21).  For  this  case  it  can  be  easily  shown 
that  the  plastic  strains  along  the  direction  defined  by  i]  ar®  zero.  In  the  an  —  an  space 
in  the  local  coordinate  axes  one  thus  obtains  strains  only  along  the  shear  direction. 

The  act\ial  incremental  plastic  strain  direction  is  obseiwed  to  ht  different  frorn  both 

•  r/--  and  n? .  To  simplify  the  issue,  one  can  assume  that  this  deviation  is  constant  in  any 

stfess  space  or  loading  level  and  measure  the  angle  of  deviation  as  6,  then  w  can  be 
measured  as  follows.  One  can  then  express  6  as  follows 


cos  (0)  =  itijr/ij 

(33) 

• 

=  +  (1  - 

(34) 

=  +  (1  “ 

(35) 

cos{9)  =  a;  -f  (1  -  w)nyny. 

(36) 

By  measuring  0  at  a  defined  loading  level  one  can  then  derive  the  value  of  u  that  would 

be  appropriate.  _  mim 

Similarly  for  loading  cases  in  the  <772 space,  one  can  choose  r\-  (0,1,0) 
(although  physically  this  does  not  represent  the  actual  fiber  direction).  This  enables  one 
to  control  the  direction  of  plastic  strains  in  this  space  and  henc»  one  can  inwrporate 
the  non-associativity  of  the  flow  rule  in  this  stress-space.  This  is  illustrated  in  Fig.  3. 


O  Yield  Surface  T 


— -  Constrained  Yield  Surface ‘g‘ 
with  n  *  (0,1 ,0) 

Fig.  3.  Illustration  of  non-associativc  flow  rule  in  022  -  012  space. 
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The  factor  u  can  be  adjusted  to  suit  experimentally  determined  values  of  the  plastic 
strain.  Constraints  in  multiple  directions  can  simultaneously  be  defined  using  one 
equation  by  appropriately  modifying  (23),  details  of  which  are  not  presented  here.  The 
plastic  strains  can  then  be  expressed  as  follows, 


<L> 

H 


(37) 


(38) 


where  L  =  a  =  Cijnij  is  the  projection  of  the  stress  rate  dy-  on  and  d  >  0  (for  plastic 
loading  after  yielding).  H  is  the  plastic  modulus.  A  detailed  description  of  the  deter¬ 
mination  of  the  plastic  modulus  for  the  anisotropic  material  is  given  later  in  this  paper. 


III.4.  Hardening  rule 

To  determine  the  shape  and  location  of  subsequent  yield  surfaces  as  loading/ 
unloading  progresses,  the  hardening  parameter  is  introduced  in  the  yield  function. 
From  the  experimental  data  of  Dvorak  et  al.  [1988]  it  has  been  demonstrated  that  the 
predominant  form  of  hardening  is  kinematic  hardening.  Also  no  significant  distortion 
of  the  shape  of  the  yield  surface  has  been  reported.  Using  these  observations  as  the 
basis  for  formulating  the  hardening  rule,  it  is  proposed  to  use  here  a  purely  kinematic 
hardening  rule  of  the  Phillips  type.  It  is  assumed  that  there  is  no  significant  change  in 
the  size  of  the  initial  yield  surface.  Hence  the  effect  of  “proportional  hardening" — 
which  accounts  for  different  expansion  of  the  yield  surface  in  different  directions — is 
ignored.  This  hardening  is  accounted  for  by  modifying  the  form  of  the  yield  surface  as 
follows. 

f—  Wu  ~  Oiij)Mijki{cia  -  ocu)  -  1.0  .  (39) 

The  evolution  equation  for  the  kinematic  hardening  rule  based  on  the  Phillips  rule  can 
be  expressed  as  follows. 


=  tl&ij  (40) 

The  above  equation  can  be  rewritten  in  the  following  form  as  follows, 

^ij  ~  (41) 

where 

lij  =  d(Tij/\\da„\\  (42) 


defines  the  unit  tensor  along  the  direction  of  loading.  The  norm  of  dy,  ||d„||,  is  found 
from  the  consistency  condition  as  follows. 


dcij  docij 


0. 


(43) 
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For  finite  increments  of  loading  and  rate  independent  plasticity  the  constitutive 
relationships  can  be  obtained  by  first  expressing  A/,  for  a  constant  M  such  that, 

A/= /(<T  +  Act,  a  +  Aa)  -/(<7,  a)  =  0  . 


Expanding  the  above  equation  one  obtains 

Mijki{(Tij  +  A(Ty  -  tty  -  Aaij){au  +  Aaw  -ocu-  Aa*/) 

-Mijkiioij  -  oii]){<Tu  -  a«)  -  1  =  0 


(44) 


Substituting  for  Aoj,  from  (41),  as  it  is  rate  independent,  it  can  be  shown  that  1|  Aa„|| 
can  be  solved  from  the  quadratic  equation 


where 


a||Aa„lP  “  i>ll  Aa„||  +  c  =  0 

b  +  <^uUj)  + 

+Mijki{aijlki  +  /ya/y) 

a  =Mijkilijlu 

c  =Mijki{ffijAau  +  Affijffu) 

+Mijki{A(TijAaki)  -  MijuiAffijau  +  aijAau). 


(45) 

(46) 

(47) 

(48) 


Once  the  magnitude  of  |1  Aa„||  is  found,  then  the  evolution  of  backstress  for  the  yield 
surface  is  obtained  and  updated  as  follows: 


^new  _  ^old  ,  a^.. 
^ij  ~  ^ij  ^  • 


(49) 


IV.  CYOJC  PLASTICITV  MODEL 

The  proposed  plasticity  model  is  further  developed  to  model  the  behavior  of  the 
composite  material  under  cyclic  loading  situations.  To  model  this  behavior,  a  modified 
form  of  the  bounding  surface  model  is  used.  The  determination  of  the  plastic  modulus 
H  is  the  main  aim  of  this  section. 

For  initially  isotropic  materials  the  plastic  modulus  is  usually  expressed  in  functional 
form  as 

jy  =  /r[l+g(6to)./(5in,«)]  (50) 

where  6i„  and  6  are  the  proximity  parameters  and  H*  is  the  asymptotic  value  of  the 
plastic  modulus.  Various  functional  forms  for  g(6J  and/(6«,  6)  have  been  proposed 
by  Dafalias  and  Popov  [1976]  and  McDowell  [1989].  Voyiadjis  and  Sivakumar  [1991, 
1994]  had  proposed  a  more  general  form  for  this  expression  as, 

/f=/r[l+g(6i„,C)./(«m,^)]  (51) 


where  the  extra  term  C  is  a  parameter  introduced  to  blend  the  Phillips’s  deviatoric  stress 
rate  direction  rule  for  the  motion  of  yield  surfaces,  when  the  yield  and  the  bounding 


surfaces  are  away  from  each  other  and  the  Tseng-Lee  [1983]  rule  for  nesting  of  yield 
surfaces,  when  the  two  surfaces  are  closer  to  each  other. 

For  isotropic  materials,  the  functional  form  of  the  plastic  modulus  as  defined  by 
Dafalias  and  Popov  [1976]  is  of  the  form 

^  =  /r+A6/(6i„-6).  (52) 

In  the  above  equation,  H*  is  the  asymptotic  value  of  the  plastic  modulus,  A  is  a  positive 
shape  parameter  and  is  the  initial  value  of  the  distance  between  the  yield  and  the 
bounding  surfaces  when  plastic  behavior  begins  in  the  loading  direction.'  For  isotropic 
materials,  the  values  of  the  three  parameters  H*,  h  and  6m  would  be  the  same,  irre¬ 
spective  of  the  location  of  the  current  stress  point  on  the  yield  surface,  and  for  loading 
along  any  direction.  This  is  not  the  case  for  anisotropic  materials. 


IV.  1.  Plastic  modulus  for  initially  anisotropic  materials 

Most  of  the  models  mentioned  above  have  been  developed  for  initially  isotropic 
materials,  although  anisotropy  is  induced  due  to  plastic  strains  subsequently.  In  such 
cases  it  is  assumed  that  both  H*  and  6m  remain  the  same  along  any  direction  (at  any 
point  on  the  yield  surface).  However  for  initially  anisotropic  or  orAotropic  materials 
the  asymptotic  value  of  the  plastic  modulus  H*  need  not,  and  in  most  cases  will  not, 
be  the  same  for  all  points  on  the  yield  surface.  For  materials  where  we  assume  the 
behavior  in  tension  and  compression  to  be  similar,  it  is  reasonable  to  assume  that  at 
mirror  image  points  of  the  yield  surface,  this  asymptotic  value  of  the  plastic  modulus  is 
the  same. 

In  Fig.  4,  which  shows  the  yield  and  the  bounding  surfaces  in  the  022  —  012  stress 
space,  points  A  and  B  are  the  location  of  the  stress  points  for  initial  yielding  for  loading 


Fig.  4.  Illustration  to  explain  plastic  modulus  determination  for  an  anisotropic  material. 
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in  the  <722  and  aij  respectively.  From  the  two  uniaxial  stress-strain  behaviors,  it  is 
observed  that  the  plastic  modulus  and  hence  the  values  of  the  three  parameters  that  are 
required  to  determine  the  plastic  modulus,  are  different. 

In  order  to  incorporate  the  above  phenomenon,  a  matrix/tensor  form  of  H  should  be 
used  to  provide  for  a  selection  of  the  hardening  parameters  corresponding  to  the  type 
of  loading  that  the  composite  is  subjected  to.  On  the  other  hand  the  form  of  (51)  that  is 
adopted  for  the  determination  of  the  plastic  modulus  is  scalar  and  has  three  parameters 
that  are  scalar  in  nature.  In  this  paper  the  authors  propose  the  following  form  for  the 
determination  of  the  parameters,  while  maintaining  their  scalar  nature. 

Consider  the  anisotropic  yield  surface  proposed  here  in  the  (erja  ~  <712)  stress 
subspace.  It  can  be  represented  as  an  ellipse  as  shown  in  Fig.  4.  One  can  assume  the 
general  location  of  the  yield  surface  as  shown.  The  location  of  a  general  stress  point  cr 
when  yielding  has  occurred — represented  by  point  C  in  the  figure  with  respect  to  the 
center  of  the  yield  surface  q,  is  represented  by  (o’  —  oc).  The  distance  of  C  from  the 
center  of  the  yield  surface  is  then  given  by 

7  =  -  av)(<^y  -  «v)  •  (53) 


A  unit  tensor  along  this  direction  can  be  written  as 

‘  (54) 

As  mentioned  earlier  the  observed  values  of  the  parameters  involved  in  the  determi¬ 
nation  of  the  plastic  modulus  are  different  along  different  loading  directions.  This 
could  be  achieved  by  using  the  tensors  in  the  form  of  second  order  tensors  Ay  and 
ITij.  These  are  then  converted  to  a  scalar  valued  form  by  taking  the  inner  product 
of  these  tensors  with  another  second  order  tensor  py  and  representing  the  result  as 
follows: 


H*  =  HlPij  (55) 

A  =  fiijpij  (56) 

6“  =  BfjPij  :  (57) 

The  expression  for  the  plastic  modulus  can  then  be  expressed  as 

•*  (58) 

where  6  is  the  distance  between  the  stress  point  on  the  yield  surface  and  the  image  point 
on  the  bounding  surface,  as  explained  later. 

Two  possible  choices  for  py  are  proposed  here  and  the  option  of  using  cither  one 
would  really  depend  on  the  observed  physical  behavior. 

IV.  1.1.  Choice  of  pij  =  Ay.  This  choice  of  py  essentially  states  that  the  values  of  the 
parameters  are  dependent  on  the  distance  of  the  location  of  the  stress  point  on  the  yield 
surface.  This  could  be  expressed  as  follows. 
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- Loading  Path 

o  Yield  Surface 

\  ,•  Bounding  Surface 


Fig.  5.  Loading  path  to  illustrate  appropriateness  of  the  choice  of  py. 


It  ensures  that  the  value  of  the  parameters  along  the  axial  directions,  assume  the 
respective  values.  At  intermediate  points  it  is  obtained  through  interpolation.  Although 
these  conditions  are  satisfied,  sometimes  this  is  not  a  very  appropriate  choice,  especially 
under  non-proportional  loading.  As  shown  in  Fig.  5,  if  the  composite  is  first  loaded 

uniaxially  along  the  cji  direction  into  the  plastic  range,  the  only  non-zero  term  of  Ay  is  • 

^22*  1.0.  This  ensures  that  the  three  values  of  the  parameters  chosen  are  those  corre¬ 
sponding  to  those  for  uniaxial  loading  along  this  direction.  If  the  material  is  then 
loaded  along  the  aji  direction,  the  values  of  the  parameters  chosen  would  be  close  to 
the  values  that  existed  for  the  previous  loading,  since  the  subsequent  location  of  the 

point  on  the  yield  surface  is  close  to  the  previous  point.  However  for  loading  along  the  • 

shear  direction  it  would  be  more  reasonable  for  the  parameters  to  assume  those  values 

corresponding  to  the  shear  direction.  Hence  this  choice,  which  has  been  tried  fiirst  in 

this  model  and  found  unsuitable  for  the  above  reason,  might  not  be  the  best  choice. 

rv.1.2.  Choice  of  Pa-hy  Since  ly  is  the  unit  tensor  along  the  loading  direction  the 

choice  of  py  =  ly  would  essentially  ensure  that  the  values  of  the  parameters  chosen  are  • 

dependent  on  the  direction  of  loading.  This  choice  is  adopted  here  in  the  proposed 
model,  since  it  has  been  found  to  be  more  suitable.  Also  this  model  is  thus  open  for 
adoption  of  any  other  suitable  choices.. 


rv.2.  Determination  of  bounding  surface 

The  initial  bounding  surface  used  here  is  an  identical  expansion  of  the  initial  yield 
surface.  This  is  done  by  determining  the  >'-intercept  of  the  bound  lines  at  the  zero 
plastic  strain  level,  in  the  stress-plastic  strain  curve.  The  bounding  surface  is  expressed 


® 


(60) 
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where  is  the  image  point  on  the  bounding  surface  and  /?(,  is  its  center  and 

'kx  0  O' 

0  JFa  0  (61) 

.0  0  £3, 

0  h  h 

bij=  h  Q  h  ‘  .  (62) 

M  h  0 . 


The  values  of  ki{i  —  1, . . . ,  6)  are  computed  from  experimental  data  along  the  respec¬ 
tive  axial  directions.  The  initial  yield  and  the  initial  bounding  surfaces  are  assumed  to 
be  congruent.  This  is  to  ensure  that  the  normals  at  the  stress  point  on  the  yield  surface 
and  the  image  point  on  the  bounding  surface  correspond  to  each  other.  For  example,  a 
uniaxial  state  of  stress  would  correspond  to  the  same  uniaxial  state  of  stress  on  the 
bounding  surface. 

The  motion  of  the  bounding  surface  is  assumed  to  be  dependent  on  the  motion  of 
yielding  surface  and  constrained  such  that  when  the  two  surfaces  intersect,  they  do  so 
tangentially.  This  is  ensured  by  finding  the  image  point  on  the  bounding  surface  from 
the  state  of  stress  on  the  yield  surface,  such  that  the  normals  to  the  two  surfaces  at  the 
respective  points  are  identical.  Given  the  state  of  stress  on  the  yield  surface  one  can  find 
the  normal  to  it  n{j  and  using  this  value  one  can  locate  the  image  point  on  the  bound¬ 
ing  surface  by  solving  a  set  of  non-linear  algebraic  equations — using  the  Newton- 
Raphson  technique.  Since  the  normals  to  the  two  surfaces  are  the  same,  n\j  =  n{j.  We 
can  express  in  terms  of  the  bounding  surface  parameters  and  image  point  stress 
values.  It  can  be  expressed  in  functional  form  as  follows: 


All  but  one  of  the  equations  represented  by  which  are  non-linear  functions  of  the 
image  point  of  stress  on  the  bounding  surface,  are  used  along  with  the  bounding  sur¬ 
face  equation  (which  is  also  a  non-linear  function  of  stresses).  All  the  equations  for  n)j 
cannot  be  used  because  they  are  interrelated  by  the  expression  n\jn\j—\. 


rV.3.  Backs  tress  for  the  bounding  surface 

The  evolution  of  the  center  of  the  bounding  surface  in  the  stress  space  as  loading 
continues  is  related  to  the  evolution  of  backstress  for  the  yield  surface  as  well  as  the 
relative  distance  between  the  stress  point  and  the  image  point.  As  given  in  Dafalias  and 
Popov  [1976]  it  can  be  expressed  as  follows: 


where 


82 


G.  Z.  Voyiadjis  and  G.  Thiagarajan 


Here  we  use  a  modified  fonn  as  follows, 


As  the  two  surfaces  touch  each  other  fiy=  0  and  hence  $ij  =  d/y  for  all  subsequent 
loadings.  From  this  increment,  the  total  backstress  for  the  bounding  stirface  is  then 

computed.  . 

The  controlling  parameter  5  in  the  determination  of  the  plastic  modulus  H  is  then 

found  from  the  formula 


V.  ELASTO-PLASnC  STIFFNESS  MATRIX 

In  order  to  determine  the  elasto-plastic  matrix  for  the  unidirectionally  oriented 
MMC  lamina,  one  has  to  incorporate  all  the  procedures  outlined  above.  Small  defor¬ 
mations  and  rate  independence  of  plastic  strains  are  assumed.  This  allows  one  to  use  an 
additive  decomposition  of  the  incremental  strain  tensor  dey  such  that, 

dey  =  d<,-hd4  (68) 

where  dc^  is  the  elastic  part  and  de^}  is  the  plastic  part  of  the  strain  tensor.  The  incre¬ 
mental  stress-strain  relations  can  be  expressed  as  follows: 

d(7y  =  Cyjfc/dcjy  (65) 

=  Cfyk/(dc/c/  ~  ficjy)  •  (^®) 

Using  (38)  for  the  plastic  strain  part  we  can  write  the  above  equation  as  follows: 

day  =  Cijkiideu  -  ^ ««)  •  (71) 

Equation  (71)  may  also  be  expressed  as  follows: 

dy  =  Ciju^ki  -  Cijkinu&/H .  (72) 

Taking  the  inner  product  of  both  sides  with  rty  one  obtains, 

CijTiij  =  Cijkinijiki  -  Cijkinijnki& /H —  a .  (73) 

Or, 

A  ^  ^  ,  (74) 


From  the  above  equation  one  can  express  d  as  follows; 
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.  Cipciltijikl  „ 
(X  ss  tl  • 


(75) 


H  +  CijuTiijnki 

Hence  the  expression  for  plastic  strains  using  (34)  can  be  written  explicitly  as  follows: 


H-\-Cpgrsn„nn 


(76) 


Substituting  this  in  the  equation  for  the  incremental  stress— strain  relations  one  obtains. 


d<Ty  *  Cijkl 


Cabcd^ab^^cd  ^ 

jr.r  n 
«  *  ^pqrs^pq^rs 


■  ^  j  ^iJklCaifcd^ab^kl^^cd 


Interchanging  the  indices  kl  with  cd  in  the  second  term  of  the  above  equation 
obtains, 

j  ^  j .  CijedCabkl^ab^a^^kl 

dCy  —  CijklC^U  ~  U  I  ^  „  m 

11  +  Lpqrs^pqnrs 


(77) 

(78) 

one 

(79) 


or, 


dir,;  = 


(80) 


where  is  the  elasto-plastic  stiffness  of  the  material  and  is  expressed  as  follows: 


Cijcd^edCgbki^ab 

H  +  Cpqrs^pqTtfs 


(81) 


VI.  NUMERICAL  SIMULATION  FOR  GENERATING  STRESS-STRAIN  CURVES 

The  above  model  has  been  built  into  a  computer  program  MC-PLAST,  for  the 
generation  of  stress-strain  curves.  Input  includes  the  elastic  constants  of  the  material, 
the  initial  yield  and  bounding  surface  parameters,  the  values  of  the  plastic  modulus 
constants  and  the  non-associativity  parameter  w.  This  program  also  reads  in  a  sequence 
of  loading  in  incremental  form  to  be  applied  and  outputs  stresses,  strains  and  plastic 
strains  at  the  end  of  each  increment. 

In  order  to  exhibit  the  validity  of  the  proposed  model,  the  model  here  is  compared 
with  results  obtained  from  the  experimental  data  of  Dvorak  et  al.  [1988]  and  Nigam  et 
al.  [1993].  The  parameters  and  constants  have  been  determined  from  these  experimental 
data.  The  experiments  had  been  done  using  boron-aluminum  tubular  composite 
specimen  having  unidirectional  lamina.  The  fibers  are  aligned  along  the  direction  of  the 
axis  of  the  tube. 

Under  plane  stress  conditions,  the  component  form  of  the  yield  surface  (7)  reduces  to 
the  following  form: 
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+  0(M2+*5)»i2-i  =  o.  ^ 

From  die  experimental  data  the  following  values  of  yield  stresses  have  been  taken. 

(rl[i  -  160.0  MPa,  =  45.0  MPa  and  =  25.0  MPa.  Using  these  data  the  following 
values  of  ^1,  kj  and  for  the  initial  yield  surface  have  been  computed  as  ki  =  0.0133, 

k2=0.0471,/k4=0.0421.  • 

A  number  of  loading  situations — both  actual  and  simulated — are  presented  here  in 
order  to  show  the  validity  and  the  stability  of  the  model. 


VI.  1.  Experimental  comparison  and  discussions 

The  loading  path  shown  in  Table  1  that  has  been  used  by  Nigam  et  al.  [1993]  has  been 
used  for  comparison.  From  the  data  they  had  presented  in  the  paper,  the  following 
values  for  the  bounding  surface  have  been  evaluated.  Cjj  =  196.0  MPa,  =  91.5 
MPa  and  =  34.0  MPa.  This  results  in  computed  values  of  the  initial  bounding 
surface  parameters  of  =  0.0108,  *2  =  0.0232  and  Ica  =  0.0323.  The  values  for  ks,  ks 
and  k^  are  not  needed  here  and  have  been  taken  to  be  zero.  However  in  order  to  make 
the  initial  yield  and  bounding  surface  congruent,  a  value  of  k\  —  0.0065  is  chosen  and  ^4  is 
modified  to  it4  =  0.0338.  This  does  not  effect  the  outcome  of  results  of  this  model  because 
this  modification  affects  the  plastic  behavior  of  the  composite  in  the  fiber  direction. 
However  for  the  MMC,  since  a  relatively  high  value  of  the  yield  stress  in  the  fiber 
direction  is  observed,  it  ensures  relatively  plastic  inextensibility  in  the  fiber  direction. 


Table  I.  Actual  experimental  loading  sequence  (and  modelled)  in 

(Nigam al.  (1993]) 


Point 

<T22 

(MPa) 

<tj2 

(MPa) 

0 

0.00 

0.00 

1 

18.00 

0.00 

2 

109.20 

0.00 

3 

1.20 

0.00 

4 

43.00 

0.00 

5 

70.00 

41.20 

6 

70.00 

33.60 

7 

1.00 

20.50 

8 

70.00 

33.60 

9 

34.10 

28.00 

10 

1.50 

-19.00 

11 

8.00 

-9.50 

12 

38.10 

-9.40 

13 

66.60 

36.00 

14 

62.00 

29.75 

15 

35.00 

29.60 

16 

6.50 

-14.43 

17 

8.00 

-9.50 

18 

38.00 

-9.50 
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The  computation  of  the  plastic  modulus  constants  is  the  next  step  in  this  process,  i.e. 
the  evaluation  of  and  H]j  and  hij.  Different  values  of  these  constants  are  evaluated  from 
experimental  results  of  the  uniaxial  stress-plastic  strain  curves  along  different  stress 


Fig.  6.  Shear  stress-plastic  shear  strain  curve  comparison  of  experimental  and  model  generated  results,  using 

non>associative  flow  rule. 


Fig.  7.  Transverse  stress-plastic  strain  curve  comparison  for  experimental  and  model  generated  results,  using 

non>associative  flow  rule. 
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directions.  The  values  of  these  constants  have  been  e\'aluated  as  follows.  /fJi  - 1,600, 
000  MPa,  =  12,000  MPa  and  /fJj^^OOO  MPa  and  the  values  of  the  other  para¬ 
meter  hij  are  An  “  9,650,000,  *22= 90,000  and  A21  =40,000  MPa,  respectively.  The  other 
values  of  this  tensor  are  assumed  to  be  zero. 


Fig.  8.  Shear  siress-plastic  shear  strain  curve  comparison  of  experimental  and  model  generated  results,  using 

associative  flow  rule. 


Fig.  9.  Transverse  stress-plastic  strain  curve  comparison  for  experimental  and  model  generated  results,  using 

associative  flow  rule. 
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In  order  to  incorporate  the  non-associativity  of  the  flow  rule  that  has  been  built  into 
the  model,  the  value  of  a>  has  been  chosen,  by  trial  and  error,  as  u;“0.5. 

Figure  7  shows  the  comparison  between  the  experimentally  obtained  and  model¬ 
generated  a22  -  curve  and  Fig.  6  shows  the  same  for  <721  -  2c^2-  From  the  com¬ 
parison  of  the  experimental  results  of  Dvorak  et  al  [1988]  and  Nigam  et  al.  [1993]  and 
the  model  generated  stress-strain  curves,  a  reasonably  good  correlation  is  observed.  In 
Figs  6  and  7  the  prediction  of  plastic  strains  for  the  case  of  uniaxial  loading  and 
unloading  in  the  <722  direction  (loading  path  1-2,  2-3)  is  good.  For  reloading  along  a 
multi-axial  path  4-5,  where  load  is  applied  simultaneously  along  the  022  and  <721 
directions,  the  model  has  been  able  to  successfully  predict  the  onset  of  yielding,  which 
in  the  authors’  opinion  is  a  very  significant  fact,  and  the  total  magnitude  of  plastic 
strains  reasonably  well. 

The  tendency  for  ratchetting  to  occur  for  cyclic  loading,  for  five  cycles  of  loading 
path  12-13-14-15-16-17  has  also  been  observed.  But  the  tendency  to  stabilize  has 
been  different  for  the  experimental  and  model  predicted  results.  This  is  because  a 
drastic  degradation  of  elastic  behavior  has  been  observed  in  the  experimental  results. 
For  the  results  of  the  experimental  behavior  the  readers  are  referred  to  the  paper  by 
Nigam  et  al  [1993].  It  has  been  unable  to  measure  this  degradation  in  elastic  behavior 
and  incorporate  them  in  the  model,  where  a  constant  elastic  behavior  is  assumed 
throughout  the  entire  loading  sequence. 

Another  important  feature  that  has  been  proposed  here  and  used  successfully  is  the 
non-associative  flow  rule  for  the  prediction  of  plastic  strains.  In  order  to  demonstrate 
that  we  must  adopt  a  non-associative  flow  rule,  the  model  is  run  with  the  same  loading 
situation,  but  with  uj  =  1 .0,  which  results  in  the  usage  of  an  associated  flow  rule.  Figures  8 
and  9  show  the  comparison  of  model  and  experimental  results  for  this  case.  For  a  pure 
associative  flow  rule  (u;=  1.0)  it  is  seen  that  plastic  strains  C22  have  been  overpredicted. 


Fig.  to.  Stress-plastic  strain  curve  for  simulated  uniaxial  cyclic  loading  of  022  stress  only. 
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A  factor  of  («;=0.5  which  incorporated  non-associativity  into  the  model  has  been 
successfully  tised  to  predict  the  plastic  strains  reasonably  in  this  direction. 

A  number  of  simulated  cyclic  loading  situations  have  been  run  on  this  model.  They 
are  uniaxial  cyclic  loading  in  the  a22  and  au  directions  and  also  radial  cyclic  loading  in 


Fig.  11.  Stress-plastic  strain  curve  for  simulated  uniaxial  cyclic  loading  for  shear  stress  only. 


# 


# 


Fig.  12.  Stress-plastic  strain  curve  for  simulated  radial  stress  loading  in  922  **  ^21- 
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the  (0-22  -  ^12)  space.  Figures  10  and  11  show  the  simulated  model  predicted  response 
under  uniaxial  cyclic  loading  in  the  <722  and  ^21  directions,  respectively.  Ratchetting 
behavior  under  constant  mean  stress  is  observed.  Figure  12  shows  the  same  for  radial 
cyclic  loading  in  the  (<722  -  <721)  stress  space.  It  is  observed  that  the  model  is  able  to 
predict  different  behavior  in  different  stress  spaces  during  the  loading  process. 

The  model  presented  here  will  be  used  in  the  micromechanical  damage  model 
proposed  by  Voyiadjis  and  Kattan  [1993].  The  proposed  model  will  be  used  to 
characterize  the  in  situ  behavior  of  the  matrix  within  the  fiber.  This  is  being  currently 
investigated  by  the  authors. 

VI.2.  Comparison  with  other  proposed  models 

The  plastic  strains  predicted  by  the  proposed  model  have  also  been  compared  with  those 
predicted  by  two  micromechanical  models,  namely  the  Periodic  Hexagonal  Array  (PHA) 
model  developed  by  Dvorak  and  Teply  [1985]  and  the  self-consistent  scheme  of  Hill  [1948, 
1965]  and  Budiansky  [1965]  using  the  Mori-Tanaka  [1973]  averaging  scheme  for  the  eva¬ 
luation  of  the  concentration  factors  by  Lagoudas  et  al.  [1991].  The  data  for  the  self-con¬ 
sistent  and  the  PHA  model  have  been  taken  from  the  paper  by  Lagoudas  et  al.  [1991]. 

This  comparison  has  been  made  in  Fig.  6,  which  shows  the  comparison  of  the  shear 
stress-plastic  shear  strain  curves  generated  by  the  models  along  with  those  from  the 
experimental  data.  It  is  seen  that  while  the  Mori-Tanaka  and  the  PHA  model  results 
underpredict  the  plastic  strains,  the  proposed  model  using  the  non-associative  flow  rule 
comes  closer  in  its  prediction.  In  the  transverse  direction  Fig.  7  the  PHA  and  the  pro¬ 
posed  model  predict  the  plastic  strains  fairly  well  whereas  the  Mori-Tanaka  model 
underpredicts  here  also.  TTie  PHA  and  the  Mori-Tanaka  model  results  have  not  been 
shown  here  to  avoid  congestion. 

vn.  CONCLUSIONS 

This  proposed  continuum  model  for  the  elasto-plastic  behavior  of  an  MMC  has  been 
able  to  successfully  capture  and  model  certain  trends  and  characteristics  that  have  been 
observed  experimentally.  Primarily  they  are,  the  usage  of  a  proposed  non-associative 
flow  rule  to  predict  plastic  strains,  incorporation  of  the  fiber-direction  plastic  inexten¬ 
sibility  criterion,  the  usage  of  a  proposed  criterion  for  evaluating  the  plastic  modulus 
and  its  incorporation  into  the  cyclic  plasticity  model.  Comparison  with  experimental 
data  has  shown  reasonably  good  correlation  and  certain  simulated  cyclic  loading 
situations  demonstrate  ratchetting  behavior.  One  feature  that  is  observed  in  the 
experimental  results,  namely  the  degradation  of  elastic  constants,  has  however  not  been 
incorporated  here  as  it  has  been  found  to  be  very  diflScult  to  measure  them. 

Adaiowkdgements—Tht  research  described  in  this  paper  is  sponsored  by  the  Air  Force  Office  of  Scientific 
Research  under  Grants  F49620-93-1-0097DEF  and  F49620-92-J0463.  The  authors  wish  to  acknowledge  the 
support  and  encouragement  of  Dr  Walter  F.  Jones.  The  authors  also  appreciate  the  experimental  data  pro¬ 
vide  to  them  by  Dr  George  J.  Dvorak  of  Rensselaer  Polytechnic  Institute.  New  York. 

REFERENCES 

1948  Hill,  R.,  “A  Theory  of  the  Yielding  and  Plastic  Flow  of  Anisotropic  Metals,"  Proc.  Roy.  Soc. 

London,  A193, 281-297. 


90 


C.  Z.  Voyiadjis  and  G.  Thiagarajan 


1951  Drucker,  D.C.,  “A  More  Fundamental  Approach  to  Plastic  Stress-Strain  Relations,”  Proc.  First 
U.S.  National  Congress  Applied  Mech.,  ASME,  New  York,  487-491. 

1954  Henhey,  A.V.,  “The  Plasticity  of  an  Isotropic  Aggregate  of  Anisotropic  Face  Centered  Cubic 
Crystals,”  J.  Appl.  Mech.  Trans.  ASME,  21,  241. 

1957  Smith,  G.  F.  and  Rivlin,  R.S.,  “The  Anisotropic  Tensors,”  25(3),  308-314. 

1959  Drucker,  D.C.,  “A  DeSnition  of  Stable  Inelastic  Material,”  J.  Appl.  Mech.  Trans.  ASME,  26, 101. 

1962  Hashin  Z.  and  Shtrikman  S.,  “A  Variational  Approach  to  the  Theory  of  the  Elastic  Behavior  of 

Polycrystals,”  J.  Mech.  Phys.  Solids,  10, 343-352. 

1965  Budiansky,  B.,  “On  the  Elastic  Moduli  of  Some  Heterogeneous  Materials,”  J.  Mech.  Phys.  Solids, 
13,223-227. 

1965  Hill,  R.,  “Continuum  micro-mechanics  of  elasto-plastic  polycrystals,”  J.  Mech.  Phys.  Solids,  13, 89- 

101. 

1967  Mulhem,  J.F.,  Rogers,  T.G.  and  Spencer,  A.J.M.,  “A  Continuum  Model  for  a  Fiber  Reinforced 
Plastic  Material,”  Proc.  Roy.  Soc.  London,  A301, 473-492. 

1969  Mulhem,  J.F.,  Rogers,  T.G.  and  Spencer,  AJ.M.,  “A  Continuum  Theory  of  an  Elastic-Plastic  Fiber 
Reinforc^  Material,"  Proc.  Roy.  Soc.  London,  7, 129-152. 

1969  Walpole,  LJ.,  “On  the  Overall  Elastic  Moduli  of  Composite  Materials,”  J.  Mech.  Phys.  Solids,  17, 
235-251. 

1970  Hutchinson,  J.W.,  “Elastic-plastic  Behavior  of  Polycrystalline  Metals  and  Composites,”  Proc.  Roy.  - 
Soc.  London,  A319, 247-272. 

1971  Huang,  W.C.,  “Plastic  Behavior  of  Some  Composite  Materials,”  J.  Composite  Mater.,  5, 320-338. 

1971  Tsai,  S.W.  and  Wu,  E.M.,  J.  Composite  Mater.,  5,  58. 

1972  Spencer,  AJ.M.,  “Deformations  of  Fiber-reinforced  Materials,”  Qarendon  Press,  Oxford. 

1973  Mori,  T.  and  Tanaka,  K.,  “Average  Stress  in  the  Matrix  and  Average  Elastic  Energy  of  Materials 
with  Misfitting  Inclusions,”  Acta  Met.,  21,  571-574. 

1976  Dafalias,  Y.F.  and  Popov,  E.P.,  “Plastic  Internal  Variables  Formalism  of  Cyclic  Plasticity,”  J.  Appl. 
Mech.,  43,  645-651. 

1979  Dvorak,  G.  J.  and  Bahei-El-Din,  YA.,  “Elastic-Plastic  Behavior  of  Fibrous  Composites,”  J.  Mech. 
Phys.  Solids,  27,  51-72. 

1983  Tseng,  N.T.  and  Lee,  G.C.,  “Simple  Plasticity  Theory  of  the  Two-Surface  Type,”  ASCE  J.  of  Engg 
Mechanics,  109,  795-810. 

1984  Eisenberg,  M.A.  and  Yen,  C.F.,  “The  Anisotropic  Deformation  of  Yield  Surfaces,”  ASME  J.  Engng 
Mater.  Technol.,  106,  355-360. 

1985  Dvorak,  GJ.  and  Teply,  J.L.,  Plasticity  Today:  Modelling,  Methods,  and  Applications,  Periodic 
Hexagonal  Array  Models  for  Plasticity  Analysis  of  Composite  Materials,  p.  623.  Elsevier,  Amster- 
dam. 

1987  Benveniste,  Y.,  “A  New  Approach  to  the  Application  of  Mori-Tanaka  Theory  in  Composite 
Materials,”  Mechanics  of  Materials,  6, 147-157. 

1987  Dvorak,  G.J.  and  Bahei-El-Din,  Y.A.,  “A  Bimodal  Plasticity  Theory  of  Fibrous  Composite 
Materials,”  Acu  Mech.,  69,  219-241. 

1988  Dvorak,  G.J.,  Bahei-El-Din,  Y.A.,  Macheret,  Y.  and  Liu  C.H.,  “An  Experiment^  Study  of  Elastic- 
Plastic  Behavior  of  a  Fibrous  Boron-Aluminum  Composite,”  J.  Mech.  Phys.  Solids,  36,  655-687. 

1988  Teply,  J.L.  and  Dvorak,  G,J.,  “Bounds  on  Overall  Instanunebus  Properties  of  Elastic-Plastic 
Composites,”  J.  Mech.  Phys.  Solids,  36,  29-58. 

1988  Dvorak,  G  J.  and  Bahei-El-Din,  Y  A.,  “Plasticity  Analysis  of  Fibrous  Composites,”  J.  Appl.  Mech., 
49,  327-335. 

1989  McDowell,  D.L.,  “Evaluation  of  Intersection  Conditions  for  Two  Surface  Plasticity  Theory,” 
Intematio^  Journal  of  Plasticity,  5, 29-50. 

1990  Voyiadjis,  G.Z.  and  Foroozesh,  M.,  “An  Anisotropic  Distortional  Yield  Model,”  ASME  J.  Appl. 
Mech.,  57,  537-547. 

1990  Weng,  G  J.  “The  Overall  Elastoplastic  Stress-Strain  Relations  of  Dual  Phase  Metals,”  J.  Mech. 
Phys.  of  Solids,  38, 3419-441 . 

1990  Zhao,  Y.H.  and  Weng,  G  J.,  “Theory  of  Plasticity  for  a  Class  of  Inclusion  and  Fiber-Reinforced 
Composites,”  Weng,  GJ.,  Taya,  M.  and  Abe,  H.  (eds).  Micromechanics  and  Inhomogeneities:  The 
Toshio  Mura  65th  Anniversary  Volume.  Springer,  New  York. 

1991  Hansen,  A.C.,  Blackketter,  D.M.  and  Walrath,  D.E.,  “An  Invariant  based  Flow  Rule  for  Ortho¬ 
tropic  Plasticity  applied  to  Composite  Materials,”  J.  Appl.  Mech.,  58(4),  881-888. 

1991  Lagoudas,  D.C.,  Gavazzi,  A.C.  and  Nigam,  H.,  “Elastoplastic  behavior  of  metal  matrix  composites 
based  on  incremental  plasticity  and  the  Mori-Tanaka  averaging  scheme,”  Comp.  Mech.,  8,  193— 
203. 

1991  Voyiadjis,  G.Z.  and  Sivakumar,  S.M.,  “A  Robust  Kinematic  Hardening  for  Cyclic  Plasticity  with 
Ratchetting  Effects.  Part  I:  Theoretical  formulation,”  Acta  Mech.,  90, 105-123. 

1992  Spencer,  AJ.M.,  “Plasticity  Theory  for  Fiber-reinforced  Composites,”  J.  Engng  Math.,  26, 107-1 18. 


# 


# 


n  wjuuv  ttujjvuvpivpiaauwij  uiv^wr 


# ' 


t 


1993  Mgam,  H.,  Dvorak,  G  J.  and  Bahei-El-Din,  YA.,  “An  Experimental  Investigation  of  Elastic-Plastic 

Behavior  of  a  Fibrous  Boron-Aluminum  Composite.  I.  Matrix-Dominated  Mode,”  Int  J.  Plasticity, 
10,23-48. 

1993  Schmidt,  RJ.,  Wang.  D.Q.  and  Hansen,  A.C.,  “Plasticity  Model  for  Transversely  Isotropic  Materials,” 

J.  Engng  Mech.,  119(4),  748-766. 

1993  Voyiadjis,  G.Z.  and  Kattan,  P.I.,  “Local  Approach  to  Damage  in  Elasto-Plastic  Metal  Matrix 
Composites,”  Int.  J.  Damage  Mech.  2(1),  92-114. 

1994  Voyiadjis,  G.Z.  and  Sivakumar,  S.M.,  “A  Robust  Kinematic  Hardemng  Rule  for  Cyclic  Plasticity 
with  Ratchetting  Effects.  Part  II:  AppUcation  to  Nonproportional  loading  Cases,”  Acta  Mech.  107, 
117-136. 

1994  Voyiadjis,  G.Z.,  Thiagarajan,  G.  and  Petralds,  E.,  “Constitutive  Modelling  for  Granular  Media 
Using  an  Arusotropic  Distortional  Yield  Model,”  Acta  Mech.,  110, 151-171. 

1995  Vo^adjis,  G.Z.  and  Thiagarajan,  G.,  “An  Anisotropic  Yield  Surface  Model  for  Directionally 
Reinforc^  Metal  Matrix  Composites,”  Int.  J.  Plasticity,  11, 867-894. 


# 


AUTHORS’  NAMES: 

George  Z.  Voyiadjis^  and  Rainer  Echle^ 

•  TITLE  OF  PAPER: 

A  Micro-Mechanical  Fatigue  Damage  Model  For  Uni-Directional  Metal  Matrix  Composites 
AUTHORS’  AFFILIATIONS: 

®  ^  Boyd  Professor,  Department  of  Civil  and  Environmental  Engineering,  Louisiana  State  University, 

Baton  Rouge  LA  70803. 

^  Doctoral  Student,  Department  of  Civil  and  Environmental  Engineering,  Louisiana  State  Univer- 

#  sity.  Baton  Rouge  LA  70803. 


ABSTRACT:  Improvement  in  design  and  enhancement  in  performance  of  aerospace  vehicles  calls 
for  the  development  of  advanced  materials  capable  of  sustaining  the  arising  loading  conditions 
while  maintaining  their  structural  integrity.  Special  consideration  has  to  be  given  to  the  behavior 
of  such  materials  under  fatigue  loading  conditions  which  are  dominating  the  flight  regime  loads.  A 
micro-mechanical  fatigue  damage  model  for  uni-directional  metal  matrix  composites  is  proposed. 
Damage  evolution  is  considered  at  the  constituent  level  through  the  application  of  the  Mori- Tanaka 
averaging  scheme.  Individual  damage  criteria  for  the  constituents  are  proposed  and  employed  to 
define  damage  evolution  equations  for  each  of  the  constituents.  Numerical  results  for  high  cycle 
fatigue  loading  are  presented  for  variations  in  material  and  model  parameters. 
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Introduction 

With  the  increase  in  performance  of  aerospace  vehicles  design  factors  such  as  weight  and  material 
strength  play  an  increased  role  in  the  design  philosophies  of  such  structures.  Along  with  such  drastic 
performance  enhancements  appropriate  light  weight  materials  need  to  be  developed  that  are  capable 
of  performing  under  such  conditions  as  those  occurring  during  flight  while  retaining  their  structural 
integrity.  Such  candidate  materials  have  been  identified  among  the  composite  materials,  especially 
in  the  area  of  metal  matrix  composites.  Special  consideration  has  been  given  to  Titanium  Matrix 
Composites  (TMC)  due  to  the  fact  that  these  materials  maintain  their  excellent  strength  to  density 
ratio  even  at  elevated  temperatures.  This  intrinsic  material  property  has  drai^Ti  attention  from 
the  turbine  engine  manufacturing  industry  for  potential  use  of  such  materials  in  advanced  zdrcraft 
turbine  engines.  The  main  reason  for  this  success  is  attributed  to  the  tremendous  reduction  in 
weight  of  the  key  engine  components  leading  to  a  possible  increase  in  engine  performance  and/or  a 
reduced  fuel  consumption.  Titanium  matrix  composites  offer  higher  mechanical  properties,  better 
dimensional  stability,  £ind  strength  retention  at  elevated  temperatures,  such  as  those  occurring  in 
turbine  engines,  as  compaired  to  their  monolithic  counterparts.  Nevertheless  the  use/employment  of 
MMC’s  and  TMC’s  still  has  major  drawbacks.  First,  the  production/manufacturing  costs  for  such 
material  are  still  high  due  to  the  special  manufacturing  processes  involved.  Second  the  employment 
of  sudi  materials  in  vital  components  of  an  aircraft  or  space  vehicle,  such  as  a  turbine  engine, 
requires  a  thorough  understanding  and  control  of  the  material  behavior  under  extreme  loading 
conditions  such  as  those  occurring  during  the  regular  service  life  of  such  vehicles.  This  calls  for  the 
development  of  material  models  which  are  capable  of  predicting  real  life  behavior  of  such  materials 
with  a  deterministic  margin  of  risk.  As  of  today  the  behavior  of  such  MMC’s  and  TMC’s  is  not  yet 
fully  understood  and  appropriate  material  models  still  lack  reliability  and  applicability  as  compared 
to  those  of  their  monolithic  counterparts.  Considerable  experimental  as  well  as  theoretical  research 
efi'ort  is  put  forward  to  better  understand  and  analytically  model  metal  matrix  composites  in  order 
to  provide  a  more  efficient  use  of  the  tremendous  potential  contained  in  these  new  materials.  In 
particular,  the  literature  lacks  a  consistent  and  systematic  approach  to  the  analysis  of  cyclic  damage 
(low-cycle  /  high-cycle  fatigue)  in  high  temperature  metal  matrix  composites.  It  is  this  area  which 
is  addressed  in  this  paper,  especially  the  consistent  development  of  a  fatigue  damage  model  for 


Tini-directional  Metal  Matrix  Composites. 


Cyclic/Fatigue  Damage  Models  in  the  Literature 

Reviewing  the  literature  on  the  subject  of  fatigue  in  engineering  materials  reveals  that  the 
explanation  of  fatigue  phenomena  and  the  prediction  of  fatigue  life  have  been  the  focus  of  im¬ 
mense  research  efforts  for  the  last  50  years.  The  two  major  analytical  approaches  used  are  the 
phenomenological  approach  and  the  crack  propagation  approach.  The  former  is  concerned  with 
lifetime  prediction  for  complex  loading  histories  using  existing  lifetime  test  data,  mostly  S  —  N  data 
for  constant  amplitude  cyclic  loading.  The  second  approach  is  concerned  with  predicting  the  growth 
of  a  dominant  crack  due  to  cyclic  load  which  is  not  the  case  for  metal  matrix  composites. 

Almost  all  of  the  known  fatigue  damage  models  for  composite  materials  are  based  on  the  models 
developed  for  their  isotropic  counterparts  (Owen  zmd  Howe  1972,  Subramanyan  1976,  Srivatsavan 
and  Subramanyan  1978,  Lemaitre  and  Plumtree  1979,  Fong  1982,  Hashin  1985,  Hwang  and  Han 
1986a,b,  Whitworth  1990).  Lack  of  theoretical  knowledge  and  sufficient  experimental  tests  on  com¬ 
posite  materials  led  to  the  application  of  known  fatigue  damage  models  to  predict  fatigue  lifetime 
of  such  materials,  despite  the  fact  that  the  fatigue  behavior  of  composite  materials  is  quite  differ¬ 
ent  from  that  of  isotropic  materials,  such  as  metals.  With  improvement  in  theoretical  knowledge 
on  composite  materials  and  experimental  equipment,  a  lot  of  studies  have  been  conducted  involv¬ 
ing  fatigue  life  and  residual  strength  degradation,  modulus  degradation  and  residual  life  theories. 
However,  it  was  soon  recognized  from  the  obtained  models  that  the  material  structure  of  such  com¬ 
posites  has  to  be  included  in  the  development  of  fatigue  damage  models  in  order  to  arrive  at  more 
feasible  and  reliable  models.  Up  to  date  there  is  no  imiversal  fatigue  damage  model  based  on  the 
micro-structure  of  the  composite  material  capable  of  predicting  fatigue  life  time  for  general  fatigue 
loading  with  reasonable  reliability. 

Arnold  and  Kruch  (1991a,b)  presented  a  phenomenological,  isothermal  transversely-isotropic 
differential  continuum  damage  mechanics  (CDM)  model  for  fatigue  of  unidirectional  composites. 
This  model  is  based  on  the  CDM  fatigue  models  for  isotropic  materials  developed  at  ONERA 
(Chaboche  and  Lesne  1988,  Chaboche  1988a,b,  1987,  Lesne  and  Savalle  1987,  Lesne  and  Cailletaud 
1987).  They  considered  the  metal  matrix  composite  as  a  pseudo  homogeneous  material  with  locallv 
definable  characteristics.  Such  local  characteristics  have  been  considered  in  the  form  of  a  directioned 


tensor  representing  the  fiber  direction.  Furthermore  the  concept  of  anisotropic  failure  surfaces 
has  been  introduced  into  the  model  based  on  deformation  theories  for  high  temperature  metal 
matrix  composites  of  Robinson  et  al.  (1987)  and  Robinson  and  Duffy  (1990).  Despite  the  rigorous 
development  the  proposed  model  has  two  major  drawbacks:  (1)  the  expensive  experimental  setup 
and  exhaustive  experiments  needed  to  obtain  the  material  parameters  used  in  the  model  equations, 
and  (2)  the  employed  scalar  measure  for  the  damage  variable.  Recently  Wilt  and  Arnold  (1994) 
presented  a  fatigue  damage  algorithm  which  employs  the  fatigue  damage  model  developed  by  Arnold 
and  Kruch  (1991a, b).  They  implemented  their  algorithm  into  the  commercial  finite  element  code 
MARC  and  used  it  to  analyze  a  cladded  MMC  ring.  Results  were  presented  on  a  qualitative  basis 
since  no  experimented  results  are  a\’ailable. 

Nicholas  (1995)  recently  reviewed  fatigue  life  time  prediction  models  for  TMC’s  which  use  fun¬ 
damentally  different  approaches.  His  investigation  showed  that  various  models  are  based  on  a  single 
parameter  and  have  limited  applicability.  Two  other  models,  a  dominant  damage  model  (Neu  1993) 
and  a  Hfe  fraction  model,  show  applicability  to  various  loading  ranges,  frequencies  and  temperature 
profiles.  Neu  (1993)  pointed  out  that  despite  the  fact  that  there  exist  several  damage  mechanisms 
it  is  possible  to  consider  the  most  dominant  ones  for  modeling  and  include  the  influence  of  others 
in  those  since  their  behavior  might  be  similar.  His  model  was  able  to  match  experimental  data  for 
isothermal  and  thermo-mechanical  fatigue  for  low-cycle  fatigue  experiments.  The  life  fraction  mod¬ 
els,  which  are  based  on  the  fact  that  fatigue  damage  accumulates  simultaneously  due  to  independent 
mechanisms,  are  able  to  model  only  specific  composite  layups  for  which  their  parameters  have  been 
calibrated.  Various  other  fatigue  investigations  have  been  performed  but  their  focus  is  on  specific 
ply-stacking  sequences  of  interest  at  the  time  of  the  investigations.  In  general  it  is  found  that  even 
though  micro-mechanical  effects  or  mechanisms  are  considered  and  incorporated  into  the  models 
there  does  not  yet  exist  a  true  micro-mechanical  fatigue  damage  model  which  considers  the  material 
behavior  and  damage  evolution  in  the  constituents  individually.  The  following  proposed  micro¬ 
mechanical  fatigue  damage  model  is  intended  to  exactly  fill  in  this  gap.  It  is  considered  as  a  first 
step  along  a  consistent  route  to  develop  a  universal  micro-mechanical  fatigue  damage  model  capable 
of  modeling  various  loading  conditions  including  thermo-mechanical  effects  as  well  as  environmental 


effects  which  occur  during  the  service  hfe  of  dynamically  loaded  composite  structures. 


Damage  Mechanics  Applied  to  Composite  Materials 

Kachanov  (1958)  pioneered  the  idea  of  damage  in  the  framework  of  continuum  mechanics.  For 
the  case  of  isotropic  damage  and  using  the  concept  of  effective  stress,  the  damage  variable  is  defined 
as  a  scalar  in  the  following  manner 


V’  = 


A-A 

A 


(1) 


where  A  is  the  effective  (net)  resisting  area  corresponding  to  the  damaged  area  A.  Using  the 
hypothesis  of  elastic  energy  equivalence  (Sidoroff  1980),  the  effective  stress  a  can  be  obtained  from 
the  above  equation  by  equating  the  force  acting  on  the  hypothetical  undamaged  area  with  the  force 
acting  on  the  actual  damaged  area. 

In  a  general  state  of  deformation  and  damage,  the  scalar  damage  variable  yj  is  replaced  by  a 
fourth-order  damage  effect  tensor  M  which  depends  on  a  second-order  damage  tensor  <t>.  In  general, 
the  effective  stress  tensor  o’  is  obtained  using  the  following  relation 


&  =  M  :  a  (2) 

where  (:)  indicates  tensor  contraction  over  two  indices.  The  nature  of  the  damage  effect  tensor  M 
is  discussed  in  the  literature  by  Voyiadjis  and  Kattan  (1992,  1993a). 

In  general  the  analysis  of  composite  materials  falls  into  two  categories.  The  first  category  consists 
of  all  approaches  that  employ  the  continuum  concept  (Talreja  1987,  Christensen  1990),  where  the 
composite  system  is  treated  as  one  continuum  and  the  equations  of  anisotropic  elasticity  are  used 
in  the  analysis.  The  second  category  encompasses  all  approaches  that  use  micro-mechanical  models 
together  with  averaging  procedures  and  homogenization  techniques  (Poursatip  et  al.  1982,  Dvorak 
and  Bahei-El-Din  1982,  1987,  Dvorak  and  Laws  1987,  Dvorak  et  aJ.  1985)  to  describe  the  material 
behavior.  In  these  models,  the  composite  is  considered  to  be  composed  of  a  number  of  individual 
phases  for  which  local  equations  are  formulated.  Employing  a  suitable  homogenization  procedure 
then  allows  one  to  analyze  the  material  behavior  of  the  entire  composite  system  based  on  the  local 
analysis. 


Dvorak  and  Bahei-El-Din  (1982,  1987)  employed  an  averaging  technique  to  analyze  the  elasto- 
plastic  behavior  of  fiber-reinforced  composites.  They  considered  elastic  fibers  with  an  elasto-plastic 
matrix.  However,  no  attempt  was  made  to  introduce  damage  in  the  constitutive  equations.  Voyiadjis 
and  Kattan  (1993b),  Voyiadjis  et  al.  (1993),  Voyiadjis  and  Kattan  (1993c)  introduced  a  consistent 
and  systematic  damage  theory  for  metal  matrix  composites  utilizing  the  micro-mechanicaJ  composite 
model  of  Dvorak  and  Bahei-El-Din  (1987).  They  introduced  two  approaches,  referred  to  in  the 
literature  as  the  overall  and  the  local  approach,  which  allow  for  a  consistent  incorporation  of  the 
damage  phenomenon  in  a  composite  material  system. 

The  overall  approach  (Kattan  and  Voyiadjis  1993)  to  damage  in  composite  materials  employs 
one  single  damage  tensor  to  refiect  all  types  of  damage  mechanisms  that  the  composite  under¬ 
goes  like  initiation,  growth  and  coalescence  of  micro-voids  and  micro-cracks.  Voyiadjis  and  Park 
(1995)  improved  the  overall  approach  by  including  and  adopting  a  general  damage  criteria  for  or¬ 
thotropic  materials  by  extending  the  formulation  of  Stumvoll  and  Swoboda  (1993)  to  MMC’s.  In 
this  improved  model  all  damage  types  are  considered  but  the  model  lacks  the  consideration  of  local 
(constituent)  as  well  as  interfaciaJ  dcunage  effects.  In  contrary  to  the  overall  approawdi  the  local 
approach  (Voyiadjis  and  Kattan  1993d)  introduces  two  independent  damage  tensors,  <f>^  and  <f)^ , 
and  hence  two  independent  damage  effect  tensors,  Af"*  and  ,  to  reflect  appropriate  damage 
mechanisms  in  the  matrix  and  fibers,  respectively.  It  is  this  latter  approach  which  is  employed  in 
the  proposed  micro-mechanical  fatigue  damage  model. 

Micro-Mechanical  Fatigue  Damage  Model 

Stress  and  Strain  Concentration  Tensors  -  In  the  derivation  of  the  model,  the  concept  of  effective 
stress  (Rabotnov  1968,  Sidoroff  1980)  is  used.  The  effective  stress  is  defined  as  the  stress  in  a 
hypothetical  state  of  deformation  that  is  free  of  damage  and  is  mechanically  equivalent  to  the 
current  state  of  deformation  and  damage.  In  a  general  state  of  deformation  and  damage,  the  effective 
Cauchy  stress  tensor  &  is  related  to  the  current  Cauchy  stress  tensor  by  the  linear  relation  given  as 
in  equation  (2).  In  the  case  of  composite  materials,  similar  constituent  (local)  stress  relations  hold 


for  the  matrix  and  fiber  stress  tensors  and  ,  respectively. 


^  (3a) 

(3b) 


where  M"*  and  are  fourth-order  local  damage  effect  tensors  for  the  matrix  and  fiber  materials, 
respectively.  The  damage  effect  tensors  M'"  and  are  dependent  on  second  order  damage  vari¬ 
ables  and  <f)^,  respectively.  These  latter  second  order  tensors  quantify  the  crack  density  in  the 
matrix  and  fibers,  respectively  (Voyiadjis  and  Venson  1995).  The  crack  density  tensors  incorporate 
both,  cracks  in  the  fiber,  matrix,  as  well  as  those  due  to  fiber  debonding.  A  complete  discussion  on 
these  tensors  is  given  in  the  work  of  Voyiadjis  and  Venson  (1995). 

In  the  proposed  model  the  matrix  is  assumed  to  be  elasto-plastic  and  the  fibers  are  assumed 
to  be  elastic,  continuous  and  aligned.  Consequently,  the  undamaged  (effective)  incremental  local 
(constituent)  constitutive  relations  are  given  by: 


da^  =  :  dl^ 


(4a) 

(4b) 


The  fourth-rank  tensors  and  are  the  undamaged  (effective)  matrix  elasto-plastic  stiffness 
tensor  and  fiber  elastic  stiffness  tensor,  respectively.  The  incremental  composite  constitutive  relation 
in  the  damaged  state  is  expressed  as  follows 


do*  =  D  :  de 


(5) 


where  de  is  the  incremental  composite  strain  tensor. 

In  order  to  arrive  at  the  local  (constituent)  relations,  given  by  equations  (4),  a  homogenization 
technique  in  the  form  of  the  Mori-Tanaka  averaging  scheme  (Chen  et  al.  1992)  is  employed.  Through 
the  use  of  the  so-called  stress  and  strain  concentration  tensors,  a  relationship  between  the  global 
applied  effective  composite  stress,  a,  and  the  local  effective  stress  in  the  constituents,  is 


# 


obtained  as  follows 

(6a) 

s=  B^  :  a 

(6b) 

where  and  represent  the  effective  stress  concentration  tensors  connecting  the  local  effective 
stresses  with  the  global  effective  stresses.  In  the  damaged  configuration  the  following  relations  are 


obtained 

3 

II 

to 

3 

(7a) 

fff  =  B^  ;  <r 

(7b) 

Combining  equations  (2),  (3),  (6)  and  (7)  one  obtains  the  relation  between  the  local  stress  concen' 
tration  tensor  and  the  local  effective  stress  concentration  tensor  as  follows 


B^  =  M-^:B^:M 


(8a) 


jpm  ^ 


(8b) 


Similar  relations  may  be  obtained  for  the  deformations  in  the  effective  (undamaged)  configuration 
as  follows 


s'"  =  A'^-.i 

(9a) 

=  A^  :e 

(9b) 

where  and  represent  the  effective  strain  concentration  tensors  connecting  the  local  effective 
strains  with  the  global  effective  strmns.  In  the  damaged  configuration  the  relations  are  given  by 


c”*  =  A^:e 

(10a) 

=  A^  :e 

(10b) 

and  furthermore 


(11a) 

A^  =  (11b) 

Effective  Volume  Fractions  -  During  the  process  of  damage  evolution  in  the  material  another 
phenomenon  has  to  be  considered.  As  damage  progresses  within  each  constituent  the  effective  load 
resisting  area/volume  changes  while  the  gross  axea/volume  remains  the  same.  Since  the  distribution 
of  forces /stresses  to  the  constituents  depends  directly  on  the  area/volume  intact  to  resist  an  applied 
force/stress  there  is  a  change  in  the  allocation  of  the  external  applied  force/stress  to  the  constituents. 
This  redistribution  of  force/stress  due  to  progressing  damage  can  be  accounted  for  by  defining  the 
so-called  effective  volume  fractions  which  are  based  on  the  updated  damage  variable  during  each 
load/stress  increment.  Expressions  for  the  effective  volume  fractions  are  given  as 

C-  = - - -  (12) 

(l_^m)  +  (l_^4)^ 


(13) 

(1  -  4>eg)  +  (1  - 

where  c[  and  c”  are  defined  as  the  volume  fractions  for  the  fiber  and  matrix  in  the  virgin  material, 
respectively.  The  expressions  for  and  4){q  are  given  as 


Hg  ~ 


4>tQ  = 


with  and  defined  as  the  critical  damage  tensors  for  the  fibers  and  the  matrix,  respectively, 
and  II  •  lli2  defined  as  the  L2  -  norm  of  the  quantity  enclosed  in  the  vertical  bars. 

Proposed  Micro-Mechanical  Fatigue  Damage  Model  -  The  proposed  fatigue  damage  criterion  g 


is  considered  as  a  function  of  the  applied  stress  ,  the  damage  parameter  <f>,  the  damage  hardening 
parameter  k,  and  a  tensor  quantity  7,  which  is  explained  below.  The  equation  for  g  is  defined  by; 

p  =  r-l  (15) 

where  J  is  defined  as 


?  =  (Xki  -  'Tki)  {Yu  -  Tit) 


(16) 


The  term  (Yki—'iki)  represents  the  translation  of  the  damage  surface  and  therefore  accounts  for  dam¬ 
age  evolution  during  cyclic  loading.  The  tensor  Y  represents  the  thermo-dynamicaJ  force  conjugate 
to  the  damage  variable  cind  is  defined  as 

(17) 

with  Cijki  =  while  the  quantity  7  can  in  principle  be  compared  to  the  backstress  in  plasticity 
theory  hence  representing  in  this  case  the  center  of  the  damage  surface  in  the  thermo-djmamical 
conjugate  force  space  Y.  Its  evolution  equation  is  given  as  follows 


y  =  - 


2  i^cd  ^abpq  ^pqkl  ^kl  d"  ^pq  ^uvpq  ^uvab  f^cd) 


Hij  —  t  (pij 


(18) 


similarly  to  the  evolution  equation  for  the  backstress  in  plasticity.  The  tensor  quantity  Wij  accounts 
for  the  anisotropic  expansion  of  the  damage  surface  and  is  given  as  follows 

v)ij  =  Uij  +  Vij  (19) 


where  the  tensor  u  is  defined  as 


Uij  =  T?(j) 


(no  sum  on  i) 


The  tensor  Vij  can  be  interpreted  physically  as  the  damage  threshold  tensor  for  the  constituent 


material  considered,  while  k  represents  the  effect  of  damage  hardening  and  is  defined  as  follows 

K=  .  d4>=  fy  :  ii>di  (21) 

7^1  JO 

Damage  hardening  is  based  on  the  increase  in  the  initial  damage  threshold  due  to  micro-hardening 
occurring  at  a  very  local  material  level  (Chow  and  Lu  1989).  The  parameter  7jj  in  equation  (16) 
adds  to  this  hardening  behavior  due  to  the  movement  of  the  damage  surfaice  in  the  direction  of 
the  evolution  f  damage.  The  remaining  variables  n,  A,-,  rjj,  material  parameters  to  be 

determined  for  each  individual  constituent.  Especially  the  form  of  the  variable  will  be  discussed 
below  in  the  numerical  implementation. 

Based  on  the  thermo-dynamical  principles  a  potential  function  for  each  constituent  is  defined 
as  (Voyiadjis  and  Kattan  1993d) 

f2  =  np  +  n‘'- Ai/- Aap  (22) 

where  11^,  /  and  g  represent  the  dissipation  energy  due  to  plasticity,  the  dissipation  energ\- 

due  to  damage,  the  plasticity  yield  surface  for  the  constituent  material  considered,  and  the  damage 
surface,  respectively.  For  loading  in  the  elastic  regime  (high  cycle  fatigue)  the  terms  involving  plastic 
dissipation  energ)'  are  neglected.  The  term  11'^  representing  the  dissipation  energy  due  to  damage 
is  given  as 


=  Yij  ^ij  +  ^k 


(23) 


Applying  the  theory  of  calculus  of  several  variables  to  solve  for  the  coefficients  Ai  and  Aj  yields 


an 

dYij 


-0 


(24) 


from  which  an  expression  for  the  damage  increment  is  obtained  as  follows 


Hence  dk2  may  be  determined  using  the  consistency  condition 


(26) 


Substitution  of  the  appropriate  terms  (equations  (18)  and  (21)  into  equation  (26)  yields 


dp  =  ^:dff+||:d</>  +  |^y:d0-c^:d«^-O 


(27) 


Replacing  d(p  with  equation  (25)  an  expression  for  dA2  is  obtained  as  follows 


Backsubstitution  of  equation  (28)  into  equation  (25)  yields  an  expression  for  the  damage  increment 
for  the  appropriate  constituent  in  terms  of  a  given  stress  increment  as 


d4>mn  ^ 


dg  dq 

doki  ^  mn 


daki 


(29) 


or 


d<pij  —  ^ ijkl 


(30) 


where 


and 


Its  —  2  ^pgkl  ^kl  +  (^pg 


dMgbcd 

d<t>rs 


(31) 


(32) 

(33) 


As  stated  elsewhere  (Stumvoll  and  Swoboda  1993)  a  damaging  state  in  a  constituent  is  given  if  for 


any  state  the  damage  criterion  is  satisfied 

S  =  0  (34) 

for  that  specific  constituent.  In  general  four  different  loading  states  are  possible 


5  <  0 

(non-damaging  loading) 

(35) 

9  =  0 

(elastic  unloading) 

(36) 

9  =  0 

(neutral  unloading) 

(37) 

9  =  0 

(loading  from  damaging  state) 

(38) 

Using  equation  (29)  the  damage  increment  per  fatigue  cycle  maybe  obtained  by  integration  over 
one  stress  cycle  as 


^4‘ij  _  f 
dN 


f^max 

min 


^ijkl 


dcki  + 


r^min 

t^max 


^ijkl  dCkl 


(39) 


where  ipijki  is  given  according  to  equation  (31).  The  dependence  of  damage  on  the  mean  stress  and 
the  amplitude  of  the  stress  cycle  is  implicitly  included  through  the  integration  of  equation  (39). 

Return  to  the  Damage  Surface 

In  the  numerical  implementation  of  the  model  it  appears  that  after  calculating  the  damage 
increment  d<j>  for  the  current  stress  increment  dtr  and  updating  all  the  appropriate  parameters 
depending  on  the  damage  variable  4>,  the  damage  surface  is  in  general  not  satisfied.  Therefore  it  is 
necessary  to  return  the  new  image  point  to  the  damage  surface  by  employing  an  appropriate  return 
criteria. 

At  the  beginning  of  the  (n  +  1)*‘  increment  we  assume  that  the  damage  surface  g  is  satisfied 


(40) 


Applying  the  stress  increment  da  (assuming  a  damage  loading)  will  result  in  a  damage  increment 
d<f>  which  will  be  used  to  update  the  values  for  k  and  7.  Checking  the  damage  surface  (equation  15) 


with  the  updated  values  for  <t,  4>,  k  and  7  will  in  general  yield 


p(n+i)  «("+!), 7!"+^)^  >  0  (41) 

where 

^(n+l)  _  ^  (42) 

(43) 

„(n+l)  _  ^(n)  ^  rf„(n+l)  (44) 

^(n+1)  _  ^(n)  ^  (f-y(n+l)  (45) 


Using  a  Taylor  series  expansion  of  order  one  expands  the  left  hand  side  of  equation  (41)  to  yield 


jl"+>) 


d<t> 


(n+l) 


dn 


(n) 


(n) 

(i«(«+i)  +  ££  d-y("+u  >  0 

01 


(46) 


Recalling  the  relationships  in  equation  (18)  and  (21)  relation  (46)  is  given  by 


p(n+i)  +  </>(") +  +  >0  (47) 

The  return  to  the  damage  surface,  hence  5'"+^  =  0,  is  now  achieved  by  adjusting  the  damage 
increment  d(f>  using  a  linear  coefficient  a  such  that 


p(n+l)  ^  +  ad0(n+l),«(n)  ^  ^y(n)  .  rf^(n+l)^yn)  +  acd<^(”+^))  =  0 

(48) 


Substitution  of  the  appropriate  expressions  for  the  derivatives  in  equation  (48)  as  well  as  equa¬ 
tion  (21)  and  (18)  and  setting  the  left  hand  side  equeJ  to  zero  allows  one  to  solve  for  the  unknown 


coefficient  a  such  that 


1 

+  H 

(ll 

(n) 

4 

y 

(49) 


Numerical  Analysis 

The  above  model  is  implemented  into  a  numerical  algorithm  and  used  to  investigate  the  fatigue 
damage  evolution  in  the  individual  constituents  of  a  uni-directionally  fiber  reinforced  metal  matrix 
composite.  No  assumption,  except  those  implicitly  included  in  the  stress  and  strain  concentration 
tensors  based  on  the  Mori-Taneika  averaging  scheme  (Chen  et  al.  1992)  are  made.  The  implemen¬ 
tation  is  performed  using  a  full  3-D  modeling  hence  avoiding  any  assumptions  to  be  made  upon 
simplification  of  fourth  order  tensors  to  two-dimensional  matrix  representation.  The  Mori-Tanaka 
averaging  scheme  is  implemented  using  the  numerical  algorithm  according  to  Lagoudas  et  aJ.  (1991). 
Only  an  elastic  analysis  is  performed  at  this  time.  Since  no  experimental  data  is  yet  available  a 
parametric  study  is  conducted  in  order  to  demonstrate  the  influence  of  various  parameters  on  the 
damage  evolution  in  the  constituents.  The  constituents  are  assumed  to  consist  of  an  isotropic  mate¬ 
rial.  The  material  used  in  the  analysis  are  given  in  Johnson  et  al.  (1990)  and  are  shown  in  Table  1. 
The  fatigue  loading  is  applied  in  the  form  of  a  sinusoidal  uni-axial  loading  given  as 


0—0  -f  CT  • 

tj  tj,meon 


where 


(50) 


^11, mean  —  bbOMPa 

and  CT  . 

ij,mean 

=  0 

OU,A  =  450  JWPo 

and  o,.^^ 

=  0 

(for  ij  ^  1) 
(for  iJ  #  1) 


For  the  numerical  integration  scheme  an  adaptive  algorithm  was  implemented  such  that  the  stress 


increments  where  taken  as 


(Jion-d^jmge  state) 

A(7y  =  IMPa  if  (damage  state) 

during  the  loading  phase  to  the  mean  stress  and 

A<To=[sm(  2^  )  sm 

(during  cyclic  loading) 

with 

IS 

II 

< 

(during  a  non-damaging  state) 

A6  = - 

900 

(during  a  damaging  state) 

for  the  cyclic  loading  phase.  Here  9  represents  simply  the  phase  angle  during  the  cyclic  loading. 
The  above  limit  values  were  adopted  based  on  a  numerical  investigation  ■which  yielded  satisfactory 
behavior  of  the  model  using  the  above  values. 

The  damage  criterion  is  evaluated  within  each  increment  and  a  return  criterion  as  described  in 
equations  (48)  and  (49)  is  applied  if  ^  10“^  Except  at  the  very  first  incident  of  damage 

this  criterion  show's  satisfactory  performance  during  the  application  of  the  return  criterion  (equa¬ 
tions  (41)  -  (49)).  The  numerical  noise  at  the  initiation  of  damage  has  been  investigated  and  it 
is  fo\md  that  a  reduction  in  the  step  size  for  the  stress  increment  reduces  the  numerical  error  ap¬ 
propriately  to  fall  within  the  specified  bounds.  This  phenomenon  is  not  observed  at  any  other 
time  during  the  analysis  (Figure  1).  It  is  attributed  to  the  point  of  discontinuity  in  the  damage 
criterion  at  the  wake  of  damage.  The  flexibility  of  the  model  is  demonstrated  through  a  parametric 
study  based  on  variations  in  the  parameters  A  and  For  the  parametric  study  the  values  of  all  the 
parameters  except  for  one  are  kept  constant  in  order  to  study  the  effect  of  a  single  parameter  on 
the  model  as  shown  in  Table  2.  The  parameters  and  account  for  the  •variation  in  the  damage 
evolution  with  respect  to  the  number  of  cycles,  especially  the  increase  in  the  damage  rate  during 
the  fatigue  life  of  a  material.  The  specific  form  of  the  parameters  and  is  obtained  from  exper¬ 
imental  curves,  such  as  those  shown  in  Figures  4  and  5,  where  the  fatigue  damage  in  the  material 


is  plotted  versus  the  number  of  applied  cycles.  Since  fatigue  damage  evolution  for  a  specific  stress 
ratio  R  is  dependent  on  the  applied  mean  stress  as  well  as  the  stress  amplitude,  such  experimental 
curves  have  to  be  obtained  for  dififerent  applied  mean  stresses  and  stress  amplitudes.  The  damage 
<l>  in  the  material  during  the  fatigue  life  maybe  obtained  by  using  the  stiffness  degradation  or  an 
equivalent  method,  such  as  sectioning  and  subsequent  SEM  evaluation  of  the  specimens  for  damage 
quantification.  Upon  inspection  of  the  obtained  experimental  curves  it  is  observed  that  basically 
three  different  regions  can  be  distinguished  during  the  fatigue  life  of  the  material  (Figures  4  and  5). 
These  different  regions  pertain  to  the  damage  initiation  phase  (Phase  I),  the  damage  propagation 
phase  (Phase  II),  and  the  failure  phase  (Phase  III).  A  distinction  for  these  regions  maybe  made  by 
specifying  bounds  in  the  form  of  the  number  of  cycles  such  as  Ni  and  N2,  as  indicated  in  Figure  4. 
This  is  done  in  general  by  visual  inspection  using  engineering  judgment  and  physical  intuition.  Us¬ 
ing  these  curves  an  evolution  equation  for  ^  with  respect  to  the  number  of  cycles  N,  the  applied 
mean  stress  Omean  and  the  stress  ratio  R  maybe  established.  For  the  current  analysis,  since  no  such 
experimental  data  are  available,  the  following  forms  for  the  parameters  and  in  terms  of  A^i 
and  N2  have  been  used  and  are  given  as 


N"  -N  m  /  N"-N\ 

(1  <N<Ni) 

(51) 

Cat  -  Co  +  ^Cl  +  ^^2 

{N"  <N<N2) 

(52) 

Co  —  N'^  J 

{N  >  N^) 

(53) 

(1  <N<N2) 

(54) 

{N>N^) 

(55) 

where 

(56) 

»  =«„+A{, 

(57) 

The  results  for  the  parametric  study  in  order  to  investigate  the  influence  of  the  model  parameter 


^0  on  the  damage  evolution  in  the  matrix  are  shown  in  Figure  3  with  aJl  other  parameters  kept 
constant.  Varying  the  value  of  the  parameter  A  and  keeping  constant  will  result  in  the  curves 
showm  in  Figure  2.  Only  the  damage  variable  <f)\\  is  shown  since  the  other  components  of  <f>  are 
equal  to  zero  or  their  value  is  smaller  by  a  magnitude  of  100.  The  reference  frames  of  the  damage 
tensor  and  the  material  system  axe  identical,  hence  “1”  representing  the  fiber  direction  while  “2”  and 
“3”  indicate  the  transverse  directions.  For  clarification  it  should  be  emphasized  that  the  plateaus 
exhibited  in  Figures  2  and  3  represent  the  unloading  phase  in  the  cyclic  loading  where  no  further 
damage  occurs. 

Two  sample  analyses  of  complete  fatigue  simulations  have  been  conducted  to  show  the  capabili¬ 
ties  of  the  developed  model.  The  result  of  such  an  analysis  for  the  damage  evolution  in  the  matrix, 
in  the  fiber  and  the  overall  composite  is  showm  in  Figure  6.  Failure  of  the  entire  composite  occurs 
due  to  fiber  failure  at  about  116000  cycles  for  the  case  of  o’u.moi  =  1000  MPo  and  a  stress  ratio 
R  =  0.1.  In  a  second  complete  fatigue  simulation  failure  occurs  at  about  217000  fatigue  cycles 
for  trii.moi  =  940  MPa  and  a  stress  ratio  R  =  0.1.  The  obtained  fatigue  life  in  the  two  cases  is 
compared  with  experimental  results  for  a  uni-directional  composite  (Johnson  1989)  as  shown  in 
Figure  7.  The  results  show  satisfactory  agreement  w’hich  establishes  the  potentiaJ  of  the  proposed 
model. 

Conclusions 

A  micro-mechanical  damage  model  for  fatigue  loading  based  on  thermo-dynamical  principles  is 
proposed.  The  model  is  applied  to  uni-directionally  reinforced  MMC’s.  Only  elastic  loading  in  the 
form  of  a  uni-axial  fatigue  loading  (in  the  fiber  direction)  is  considered  hence  reflecting  high  cycle 
fatigue  loading.  Numerical  results  from  the  parametric  study  show  the  influence  of  various  model 
parameters  on  the  damage  evolution  in  the  constituents.  A  sample  analysis  for  a  complete  fatigue 
simidation  with  final  failure  is  shown. 
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Table  1:  Material  properties  used  in  the  analysis 


E  {GPa) 

1/ 

Ou  {MPa) 

oy  (MPa) 

c  (in  %) 

Matrix  {Ti  —  15  —  3) 

92.4 

0.35 

933.6 

689.5 

67.5 

Fiber  {SCS  -  6) 

400.0 

0.25 

N/A 

N/A 

32.5 

Table  2;  Model  parameters  used  in  the  analysis 


V  {MPa) 

A  {MPa) 

V 

c{MPa) 

Matrix  {Ti  —  15  —  3) 

0.1 

80000 

1.0 

refer  to  Eqs.  (51)  -  (57) 

1.0 

1.0 

6 

Fiber  {SCS  -  6) 

3 

160000 

1.0 

refer  to  Eqs.  (51)  -  (57) 

1.0 

1.0 

6 

Ni 

N2 

^0 

6 

6 

Figure 

Matrix  {Ti  —  15  —  3) 

10 

110000 

0.55 

0.02 

0.03 

6 

Fiber  {SCS  -  6) 

N/A 

110000 

0.56 

N/A 

0.03 

6 
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Validation  of  Employed  Return  Criterion  for  Damage  Surface  g 


Figure  1:  Validation  of  the  employed  return  criterion 


Influence  of  k  on  damage  evolution 


Figure  2:  Variation  in  damage  evolution  for  various  values  of  X 
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Figure  3:  Variation  in  damage  evolution  for  various  values  of 


R  =  const. 
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Figure  4:  4>  —  N  diagrams  for  determination  of  ^  for  constant  R 


Om  —  const. 

R\  >  R2  >  Rz  ^  -^4 


Figure  h:  (f>  —  N  diagrams  for  determination  of  i  for  constant  Omean 


Damage  Evolution  During  HCF  Simulation 


Figure  6:  Fatigue  damage  evolution  during  a  complete  simulation 
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Abstract— A  microroechanical  damage  composite  model  is  used  here  such  Aat  separate  1<^1  ev^tion 
damage  relations  are  used  for  each  of  the  matrix  and  fiber.  In  addition,  ^s  is  coupled  with  interfacial 
damage  between  the  matrix  and  fiber  exclusively.  An  overall  response  is  linked  to  these  damage  re¬ 
lations  through  a  certain  homogenization  procedure.  A  finite  element  analysis  is  used  for  quantifjn^ 
each  type  of  damage  and  predicting  the  failure  loads  of  dog-bone  shaped  specimen  and  center-cracked 
laminate  meul  matrix  composite  plates.  The  development  of  damage  zones  and  the 
sponse  are  shown  for  two  types  of  laminated  layups,  a  (0/90)s  layup  and  a  (±45)$  layup.  ©  1997  tlse- 
vier  Science  Ltd.  All  rights  reserved. 


INTRODUCTION 

Damage  and  plastic  deformation  is  incorporated  in  the  proposed  model  that  is  used  for  the 
analysis  of  fiber-reinforced  metal  matrix  composite  materials.  The  proposed  micromechanical 
damage  composite  model  used  here  is  such  that  separate  local  constitutive  damage  relations  are 
used  for  each  of  the  matrix  and  fiber.  This  is  coupled  with  the  interfacial  damage  between  the 
matrix  and  fiber  exclusively.  The  damage  relations  are  linked  to  the  overall  response  through  a 
certain  homogenization  procedure.  Three  fourth-order,  damage  tensors  Af"*,  m  and  Aj  are 
used  here  for  the  two  constituents  (matrix  and  fibers)  of  the  composite  system.  The  matrix 
damage  effect  tensor  A/™  is  assumed  to  reflect  all  types  of  damage  that  the  matrix  material 
undergoes  such  as  nucleation  and  coalescence  of  voids  and  microcracks.  The  fiber  damage  effect 
tensor  is  considered  to  reflect  all  types  of  fiber  damage  such  as  fracture  of  fibers.  An  ad¬ 
ditional  tensor  is  incorporated  in  the  overall  formulation  that  represents  interfacial  damage 
between  the  matrix  and  fiber.  An  overall  damage  effect  tensor,  A/,  is  introduced,  that  accounts 
for  all  these  separate  damage  tensors  Af™,  A/  and  . 


THEORETICAL  PRELIMINARIES 

The  metal  matrix  composite  used  in  this  work  consists  of  an  elasto-plastic  ductile  metal 
matrix  reinforced  by  elastic  aligned  continuous  fibers.  The  composite  system  is  restricted  to 
small  deformations  with  infinitesimal  strains.  In  the  initial  configuration,  Co,  the  composite  ma¬ 
terial  is  assumed  to  be  undeformed  and  undamaged.  The  initial  matrix  and  fiber  subconfigura¬ 
tions  are  denoted  by  C?  and  Ct,  respectively.  Due  to  applied  loads,  the  composite  material  is 
assumed  to  undergo  elasto-plastic  deformation  and  damage,  and  the  resulting  overall  configur¬ 
ation  is  denoted  by  C.  The  resulting  matrix  and  fiber  local  subconfigurations  are  denoted  by  C*" 
and  d,  respectively.  Damage  is  quantified  using  the  concept  proposed  by  Kachanov[l]  whereby 
two  kinds  of  fictitious  configurations  C  and  C  of  the  composite  system  are  considered.  C  con¬ 
figuration  is  obtained  from  C  by  removing  all  the  damages,  while  C  configuration  is  obtained 
from  C  by  removing  only  the  jnterfacial  damage  between  the  matrix  and  fiber.  C  is  termed  full 
effective  configuration,  while  C  the  partial  effective  configuration. 

A  coupling  formulation  of  plastic  flow  and  damage  propagation  seems  to  be  impossible, 
due  to  the  presence  of  the  two  different  dissipative  mechanisms  that  influence  each  other.  For 
example,  the  position  of  slip  planes  affects  the  orientation  of  nucleated  microcracks.  One  can, 
however,  assume  that  the  energy  dissipated  in  the  yielding  and  damaging  processes  is  indepen¬ 
dent  of  each  other  and  apply  a  phenomenological  model  of  interaction.  Use  will  be  made  of  the 
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concept  of  effective  stress  [2].  Assuming  a  fictitious  undamaged  system,  the  dissipation  energy 
due  to  plastic  flow  in  this  undamaged  system  is  assumed  to  be  equal  to  the  dissipation  energy 
due  to  plastic  flow  in  the  real  damaged  system. 

The  basic  feature  of  the  approach  presented  here  is  that  local  effects  of  damages  are  con¬ 
sidered  whereby  these  effects  are  described  separately  by  the  matrix,  fiber  and  interfacial 
damage.  It  is  clear  the  local  nature  of  damage  of  this  approach  whereby  the  different  damages 
are  separately  isolated.  This  approach  can  be  summarized  in  the  following  three  steps^  First, 
apply  the  local  damage  effect  tensors  and  to  the  local  effective  configurations  C“  and 
O,  respectively.  This  is  followed  by  applyinj  the  damage  stress  concentration  factors  J"  and 
to  the  local  partial  effective  configurations  C™  and  in  order  to  obtain  the  overall  partial  effec¬ 
tive  configuration  C.  Finally,  one  applies  the  interfacial  damage  effect  tensor  to  the  overall 
partial  effective  configuration  C  to  obtain  the  overall  damaged  configuration  C. 

In  the  formulation  of  this  work,  quantities  based  on  the  full  effective  configuration,  C,  and 
the  partial  effective  configuration,  C,  are  denoted  by  a  superposed  bar  and  superposed  tilde,  re¬ 
spectively.  Fiber  and  matrix  related  quantities  are  denoted  by  a  superscript  m  or  f,  appropri¬ 
ately.  In  addition,  interfacial  damage  related  quantities  are  denoted  by  a  superscript  d. 

In  the  undamaged  effective  configuration  C,  the  overall  effective  stress  increment  da  can  be 
expressed  in  terms  of  the  phase  effective  stress  increments  d5®  and  da^  as  follows: 

ddy^?^da^  +  ^dd^j  (1) 

where  c®  and  ?  are  the  matrix  and  fiber  volume  fractions,  respectively.  The  local-overall  re¬ 
lations  for  the  stress  increments  for  the  matrix  and  fibers  in  the  fictitious  local  and  overall  con¬ 
figurations  are  given  as  follows: 

da'j  =  da*/*  where  r  =  m*  f *  (2) 

where  is  the  stress  concentration  tensor  for  the  matrix  or  fibers.  Similarly,  the  overall  effective 
strain  increment  is  assumed  in  the  effective  configurations  such  that 

d4i,  =  c“d€”-hcfd4,  (3) 

where  dl“  and  df  are  the  effective  matrix  and  fiber  strain  increments,  respective!}^.  The  additive 
decomposition  is  assumed  of  the  matrix  and  overall  strain  increments  in  C°‘  and  C,  respectively, 
such  that 


dlj^diy'  +  dlj"  (4) 

d€^  =  dl?'  +  d€^^  (5) 

where  '  indicates  the  elastic  and  *  indicates  the  plastic  part  of  the  increment.  Equations  (4)  and 
(S)  are  justified  in  view  of  the  assumption  of  small  strains.  The  fibers  are  assumed  to  deform 
elastically  and  therefore  we  have 

d4  =  d4'.  (6) 

The  local-overall  relations  for  the  effective  strain  increments  are  given  as  follows: 

d^y  =  A'yutu,  where  r  =  m,  f,  (7) 

where  A'  is  the  strain  concentration  tensor  for  the  matrix  or  fibers. 

The  expressions  for  the  stress  and  strain  tensors  based  on  the  Mori-Tanaka  method  are 
given  in  Appendix  A. 

CONSTITUTIVE  EQUATIONS  FOR  THE  UNDAMAGED  METAL  MATRIX  COMPOSITE 
The  constitutive  relations  for  the  matrix  and  fibers  are  given  by  the  following  relations: 
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where  ^  is  the  elasto-plastic  stiffness  tensor  of  the  matrix  and  if  is  the  elastic  stiffness  tensor 
of  the  fiber.  The  elastic  stiffnesses  of  the  matrix  and  fiber  are  given  by  the  following  expressions 

=  (10) 

=  (11) 

where  iF,  Cm,  if  and  &  are  Lame’s  constants  for  the  matrix  and  fibers,  respectively. 
Substituting  for  da°‘  and  d^  from  eqs  (8)  and  (9),  respectively,  into  eq.  (1)  and  making  use  of 
eq.  (7),  one  obtains  the  relation 

ioy  ®  iiypq  (12) 

where 

D  is  the  elasto-plastic  stiffness  of  the  composite  in  the  effective  undamaged  configuration. 

PLASTIC  BEHAVIOR  OF  THE  UNDAMAGED  METAL  MATRIX  COMPOSITE 

The  elasto-plastic  constitutive  model  for  the  matrix  is  based  on  the  von  Mises  type  yield 
function a“)  in  the  local  configuration  C®  such  that 

r  m|(c“-a“)(5?-fij)-5f’s=0.  (14) 

^  is  a  material  constant  denoting  the  uniaxial  yield  stress  of  the  matrix  material  and  d®  is  the 
backstress  tensor.  The  plastic  flow  in  the  configuration  C™  is  given  by  the  associated  flow  rule 
in  the  form 


=  (15) 

where  A"  is  a  scalar  function  introduced  as  a  Lagrangian  multiplier  in  the  constrain  thermodyn¬ 
amic  equations  for  the  matrix  material. 

In  order  to  describe  kinematic  hardening  for  the  matrix,  the  Prager-Ziegler  evolution  law 
[3]  is  used  here  in  the  configuration  C®,  as  follows: 

d«^  =  M"(5?-fi^).  (16) 

where  pF  is  a  scalar  function  given  in  terms  of  Am  as  follows; 


8/"  3/° 


and  is  a  material  parameter. 

The  parameter  A*"  is  obtained  from  the  consistency 

d/"(5“,d“)  =  0 

and  is  given  by  the  following  relation: 

jm  1  ™ 

^  9“  BoF 


(17) 


(18) 


(19) 
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where 


a/“  3/« 

Iff:. 


The  resulting  elasto-plastic  stiffness  tensor  for  the  matrix  material  in  the  effective  undamaged 
configuration  is  given  by 

Am  _  fm  ^  m  pm  r2n 


Damage  effect  tensors 

_The  effective  stress  concept  [1,4]  is  used  in  this  work.  Considering  the  overall  configurations 
C,  C  and  C,  one  can  introduce  an  overall  damage  effect  tensor  M  and  a  partial  damage  effect 
tensor  M  for  the  whole  composite  system.  These  tensors  are  defined  similarly  to  the  definitions 
of  A/“,  hf  and  such  that 

ay  =  MyklOkl  (22) 


resent 


The? 
cedur 
by  Vc 


Oy  =  Mykiakl-  (23) 

The  tensor  M  reflects  all  types  of  damage  that  the  composite  undergoes  including  the  damage 
due  to  the  interaction  between  the  matrix  and  fibers  while  the  tensor  M  reflects  damage  of 
matrix  and  fibers,  excluding  the  interfacial  damage.  A  matrix  representation  was  expUcitly  de¬ 
rived  for  this  fourth  order  tensor  by  expressing  the  stresses  in  vector  form.  The  tensor  M  was 
shown  to  be  symmetric.  The  symmetry  property  of  the  tensor  M  is  used  extensively  in  the  deri¬ 
vation  that  follows.  The  same  holds  true  for  the  tensors  3/^  and  Similar  to  tensor 
both  tensors  and  Af  could  be  represented  in  terms  of  second  order  tensors  and  re¬ 
spectively.  The  effective  matrix  stress  and  the  corresponding  fiber  stress  are  defined  as  follows: 

=  (24) 


4  =  (25) 

where  5“  and.o*^  are  the  partial  effective  stresses  in  the  and  configurations,  resi>ectivcly. 
These  stresses  are  termed  partial  effective  since  the  interfacial  damage  has  not  yet  been  incorpor¬ 
ated  into  the.  formulation.  The  overall  damage  effect  tensor,  M,  can  be  related  to  the  partial 
damage  effect  tensor,  M,  and  the  local  damage  effect  tensors  such  as 


where 


Miikl  —  «i 


Myr,  =  ?"3/  „  -h  L5 


iipy  prr’ 


This  expression  defines  the  cumulative  damage  of  the  composite  as  a  function  of  the  local 
matrix  and  fiber  damages  Af™  and  hf,  respectively,  as  well  as  the  interfacial  damage  3f^[5,6]. 

The  expression  of  the  fourth  order  local  damage  effect  tensor  A/“,  and  can  be  rep¬ 
resented  by  a  6  X  6  matrix  as  a  function  of  (^(/  —  ^^)  in  the  form 


IMr  =  [3/(ij,  -^^)].  r  =  m,  f.  d. 


to  n# 


comp 

Thed 


where  Su  is  the  Kronecker  delta. 


T 

from 
The  I 
singit 
above 


(28) 
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In  conjunction  with  the  matrix  form  of  M  given  by  eq.  (28),  the  local  stress  tensor  is  rep¬ 
resented  by  a  vector  given  by 

[al  ala'll,  022,  =  (29) 

The  symmetrized  o'  and  o  used  here  are  given  by[4] 

Oy  =  where  r  «  m.  f  (30) 

and 


Oy  =  +  (^11  “  4il)~^Oy]. 


(31) 


The  stresses  given  by  eqs  (30)  and  (31)  are  frame-independent.  Using  the  symmetrization  pro¬ 
cedure  outlined  by  eqs  (30)  and  (31),  the  corresponding  6x6  matrix  form  of  tensor  M  is  given 
by  Voyiadjis  and  Kattanl?]  as  follows: 


2(022(033  “  2^3 
0 
0 

^13^3  +  ^12tt>33 


0 


4>12^3  +  ^13(022 


0 

2<Ull£t>33  -  2^3 
0 

^13d*23  +  ^12‘V33 
4>12^\3  +  4^23(0  U 
0 


0 

0 

2(0i\(022  -  24^2 

-P. 

^12^13  +  d»23<*>ll 
^12^3,+ ^13^^22 


2^13^3  +  2^12£i^33 
2d’I3^3  +  2^12tt>33 
0 

<022(033  +  (0U(033  —4^3~  ^13 
^12^3  +  ^13i‘>22 
^12013  +  ^3^^11 


0 

2012013  +2023tVii 
2012013  +  20230)11 
012023  +0130)22 
0)110)33-1-0)110)22  —  013  —  012 
013023  +0120)33 


2012023  +  20130)22 
0 

2012023  +  20130)22 
012013  +  0230)11 
013023  +  0120)33 
0)220)33  +  0)|iO)22  —  023  ~  012 


r  =  m,  f ,  d 


(32) 


and  V  is  given  by 

V  ss  0)110)220)33  ~  0230)11  “  0130)22  ~  0120)33  "•  2012023013-  (33) 

(Oy  is  used  to  denote  Sy-4>ij-  4>y  used  in  eq*  (32)  and  (33)  represents  0",  4>y  or  4>t  with  respect 
to  matrix  damage,  fiber  damage  or  interfacial  damage,  accordingly. 


Damage  stress  and  strain  concentration  tensors 

The  matrix  and  fiber  stress  concentration  factors  are  defined  as  fourth-rank  tensors.  As 
composites  undergo  damage  the  stress  and  strain  concentration  factors  do  not  remain  constant. 
The  damage  stress  and  strain  concentration  tensors  are  given  by  the  following  relations[6]: 

=  M  ylBr^rsMrskt,  T  =  m,  f  (34) 


^yki  ”  ^  ypi^pfrs^'Oci,  r  =  m,  f. 


(35) 


ANISOTROPIC  DAMAGE  ANALYSIS 

The  damage  mechanism  for  each  of  the  constituents  of  the  composite  material  is  different 
from  the  other.  The  matrix  undergoes  ductile  damage  while  the  fiber  undergoes  brittle  damage. 
The  mechanism  of  interfacial  damage  is  dependent  on  the  fiber  direction.  It  is  clear  that  one 
single  damage  micro-mechanism  cannot  be  considered  for  the  three  types  of  damages  outlined 
above.  We  therefore  consider  each  damage  evolution  separately. 
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Damage  criterion 

An  anisotropic  damage  criterion  is  proposed  in  this  work.  In  order  to  obtain  a  damage  cri¬ 
terion  for  nonproportional  loading,  the  anisotropy  of  damage  increase  gardening)  must  be  con¬ 
sidered.  This  is  accomplished  by  expressing  the  damage  criterion  in  terms  of  a  tensonal 
hardening  parameter,  h.  The  damage  criterion  used  here  is  of  the  form  suggested  by  Mro2[8] 
such  that 

^  E  l'(y,  A)  =  0,  r  =  m,  f,  d,  (36) 

where  is  a  generalized  thermodynamic  force  conjugate  to  the  damage  tensor  ^  for  Mch  of 
the  damages  associated  with  the  matrix,  fiber  and  debonding.  Equation  (36)  is  an  isotropic  func¬ 
tion  of  tensors  Y  and  h'  such  that 


where 


and  p'ijkt  is  equivalent  to  Hill’s  tensor  for  yield  surfaces.  The  hardening  tensor  h  is  given  by 
Tensors  i/  and  are  here  defined  for  orthotropic  materials  as  follows: 


mecb 

energ 

pov# 


The^ 


where  r  =  m,  f,  d. 


These  tensors  x/  and  Y  are  generalizations  to  orthotropic  materials  of  the  scalar  forms  for  iso¬ 
tropic  materials  originally  proposed  by  Stumvoll  and  Swoboda[9].  In  expressions  (40)  and  (41), 
the  scalar  quantities  Xj,  XJ,  Xj,  v^,  vj,  {i,  {5,  ^j,  ijt,  ri2  and  ^5  material  parameters 
obtained  by  matching  the  theory  with  experimental  results.  The  parameters  X'’’,  X'’,  X'*,  vj,  vj 
and  v5  are  explicitly  related  to  the  physical  properties  of  the  materiaU6]. 

In  eq.  (41),  V2  and  vj  define  the  initial  threshold  against  damage  for  the  orthotropic  ma¬ 
terial.  It  is  obtained  from  the  constraint  that  the  onset  of  damage  conesponds  to  the  stress  level 
at  which  the  virgin  material  starts  exhibiting  nonlinearity. 

Referring  to  eq.  (40), is  a  scalar  hardening  parameter  given  by 


=  ^^Yyd<l>'y, 


where  r  =  m,  f,  d. 


As  outlined  by  Stumvoll  and  Swoboda^],  the  damaging  sute  is  any  sute  that  satisfies  g  -  0. 
Four  states  are  outlined  here 

g^<0  (elastic-unloading)  (43) 


Maid 


g'  =  0,  dY'y<0  (elastic-unloading) 
oY  g 
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l'  *  0,  dy  «  =  0  (neutral  loading) 
or  I  ^ 


(45) 


f^  =  0, 


(loading  from  a  damaging  state). 


(46) 


In  this  section,  the  anisotropic  damage  criterion  g  is  deSned  by  eq.  (37)  as  well  as  the  loading 
conditions  outlined  by  eqs  (43)-(46).  The  anisotropic  damage  criterion  is  defined  through  the 
second  order  tensors  i/  and  V',  and  the  damage  tensor  4^'  for  each  constituent  of  the  composite 
material.  In  this  work,  we  assume  that  the  matrix  and  the  fiber  are  isotropic  materials  while  the 
anisotropic  damage  criterion  is  used  to  describe  the  interfacial  damage. 


Damage  evolution  of  the  matrix 

The  metal  matrix  exhibits  two  energy  dissipative  behaviors.  Although  the  two  dissipative 
mechanisms  of  plasticity  and  damage  influence  each  other,  in  this  work,  it  is  assumed  that  the 
energy  dissipated  due  to  plasticity  and  that  due  to  damage  are  independent  of  each  other.  The 
power  of  dissipation  for  the  matrix  is  given  by 

fi"  =  ft”^  +  (47) 

where  ft  ^  is  the  plastic  dissipation  and  n™**  the  corresponding  damage  dissipation.  The  plastic 
dissipation  is  given  by 

=  of  dc^".  (48) 

In  this  work,  a  small  strain  theory  is  assumed  and  the  strain  rate  is  assumed  to  be  decomposed 
into  an  elastic  component  and  a  plastic  component  c®",  such  that 

df7  =  d^  +  dl^".  (49) 

The  associated  damage  dissipation  is  given  by 

n"‘’  =  y?d^”  (50) 

where  T”  is  a  generalized  thermodynamic  force  conjugate  to  the  damage  tensor  0”.  The  ficti¬ 
tious  undamaged  material  is  characterized  by  the  effective  stress  and  effective  strain.  Since  in  the 
full  effective  configuration,  C™,  the  matrix  has  deformed  with  no  additional  damage,  the  dissipa¬ 
tion  energy  in  C™  is  only  composed  of  the  plastic  dissipation 

n”  =  n“’’  (51) 

and  therefore 


n"‘  =  55‘d7 


(52) 


This  is  ^cause  plastic  yielding  is  assumed  to  be  independent  of  the  damage  process.  The  plastic 
dissipation  in  the  damaged  matrix  is  equal  to  the  corresponding  plastic  dissipation  in  the  full 
effective  configuration,  C®.  This  concludes  that 


.n.p^.n>p^ 

which  implies  that 

a?‘d7''  =  o^d€“'. 

Making  use  of  eq.  (54)  together  with 

one  obtains  a  transformation  equation  for  the  plastic  strain  rates  such  that 


(53) 

(54) 


(55) 
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Making  use  of  the  calculus  of  functions  of  several  variables,  one  introduces  two  Lagrange 
multipliers  A”  and  A”  in  order  to  form  the  function  SI”'  such  that 

n“  =  n“  -  Aff”'  -  A5‘g®. 

In  eq.  (57),/"(<5®,ct®)  is  the  plastic  yield  function  of  the  matrix  and  a®  is  the  backless  tensor. 
I®  is  the  damage  potential  which  is  a  function  of  y“.  To  extremize  the  function  Cr^,  one  uses 
the  necessary  conditions 

(58) 

daf 


which  give  the  corresponding  plastic  strain  rate  and  damage  rate  evolution  equations,  respect¬ 
ively: 

=  (60) 

°  ‘J  '  35? 


A/fP  —  A“ 

UJ  y 

Equation  (61)  gives  the  increment  of  damage  from  the  damage  potential  g®.  Using  the  consist¬ 
ency  condition  for  the  matrix  damage  g® 

dg®  =  0  (62) 

one  obtains  the  parameter  Af.  Equation  (62)  states  that  after  an  increment  of  damage,  the 
volume  element  again  must  be  in  a  damaging  state.  From  eq.  (62),  one  obtains 

j  y  m 
aym®^  1) 

(63) 

^2  -  3g“  3g“  • 


34^BYf 


Substituting  eq.  (63)  into  eq.  (61),  one  obtains 


where  iji'”  is  a  fourth  order  tensor  defined  as 

3g®  3g“ 

JYf^  (65) 

3^3y» 

The  generalized  thermodynamic  force  T®  is  assumed  to  be  a  function  of  the  elastic  component 
of  the  strain  tensor  e®'  and  the  damage  tensor  (#>®,  or  the  stress  5®  and 

y»  =  y®(fj',^por  y^  =  y®(55'.«?).  (66) 

The  evolution  equation  for  1®  may  be  expressed  as  follows: 


Dar^ 

1 

and  t 


and# 

EFM  56/4 
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(56) 

Lagrange 

^  (57) 

ess  tensor. 
.  one  uses 


.  (59) 
4,  respert- 

(60) 


o 


(61) 


It  consist- 

(62) 
-sage,  the 


(63) 


fly®  fly® 


(67) 


Substituting  for  dy“  from  eq.  (67)  into  eq.  (64),  one  obtains  the  evolution  equation  for  such 
that 


(68) 


or 


where 


7“ 

^yki 


1  dy 

=  2(^ikSji  +  Mj*)  - 


(69) 


(70) 


The  thermodynamic  force  associated  with  damage  is  obtained  using  the  enthalpy  of  the 
damaged  matrix  where 

=  ^  (71) 

In  eq.  (71),  £“  is  the  damaged  elastic  stiffness  of  the  matrix.  The  thermodynamic  force  of  the 
matrix  is  given  by 


(72) 


Using  the  energy  equivalence  principle  [10],  one  obtains  a  relation  between  the  damaged 
elastic  compliance,  £”®,  for  the  matrix  and  its  corresponding  undamaged  elastic  compbance 
such  that  [11] 

£^7(^")  =  A/“i;(^”)£-“MS*;(<^").  (73) 

Making  use  of  eqs  (71)  and  (72),  the  thermodynamic  force  for  the  matrix  is  obtained  expbcitly 
such  that 

1  fl  A/  ® 


#  m 


(65) 

omponent 


(66) 

# 


Damage  evolution  of  the  fiber 

The  gradual  degradation  of  the  elastic  stiffness  of  the  fiber  is  caused  only  through  damage 
and  therefore  no  plastic  dissipation  occurs.  We  therefore  have 


n'  =  n"  =  y[.d^- 

(75) 

and 

n^  =  o. 

(76) 

Accordingly,  the  function  is  given  by 

> 

< 

1 

•C 

II 

(77) 

and 

EFMM/4-B 
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d^f.  s 

flVf  • 


Using  the  consistency  condition  for  the  damage  of  the  fiber 

d/  =  0  (79) 

one  obtains  the  evolution  equation  for  <l>* 

^  yU^kl’  (80) 

where  ^  is  a  fourth  order  tensor  similar  to  A®  expressed  by  eq.  (69).  is  obtained  in  a  similar 
approach  to  that  of  the  matrix,  y™ ,  and  has  a  similar  form,  except  replacing  the  superscript  m 
with  f. 

Jnterfacial  damage  evolution 

The  interfacial  damage  can  be  defined  in  terms  of  a  second  order  symmetric  tensor  such 
as[6] 

^^*^“(5.5).  (81) 

More  elaborate  interfacial  damage  expressions  could  be  derived  based  on  the  work  of  Levy[12]. 

The  corresponding  power  of  dissipation  due  to  interfacial  damage  is  given  by 

n^  =  Yld<pl  (82) 


n‘'  =  o. 


The  function  n"*  is  expressed  as 


S'*  =  n**  -  A  V 


d(^‘’  =  A** 

ayd- 


Using  the  consistency  condition  for  the  interfacial  damage 

dg^  =  0  (86) 

one  obtains  the  evolution  expression  for  (p^  such  that 

<i4>y  =  Xlu^ki.  (87) 

Similar  to  the  procedure  outlined  for  the  other  two  types  of  damages  y**  could  be  easily 
obtained  accordingly,  such  as 


CONSTITUTIVE  MODEL  FOR  THE  DAMAGED  MATERIAL 

Derivation  of  the  elasto-plastic  constitutive  model  for  the  damaged  composite  system  is  per¬ 
formed  in  three  steps.  The  first  step  involves  the  derivation  of  separate  constitutive  equations 
for  the  matrix  and  fiber  in  their  respective  damaged  configurations  and  C^.  This  is  followed 
by  the  second  step  which  combines  the  two  constitutive  equations  into  one  for  the  overall  com¬ 
posite  system  in  its  partial  effective  configuration  C.  FinaUy,  the  interfacial  damage  is  incorpor- 
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(78) 

#  (79) 

(80) 

1  a  similar 
^rschpt  m 

It  such 

^  (81) 
LevyI12]. 

(82) 

# 


ated  is  order  to  obtain  the  final  constitutive  equation  that  includes  all  three  types  of  damage  in 
the  damaged  configuration  C.  _  _ 

The  local  damaged  elastic  stiffness  tensors  £“  and  Ef  in  the  subconfigurations  C®  and 
respectively,  are  obtained  using  the  hypothesis  of  complemental  energy  equivalence  such 

that 

(89) 

and 

= 


In  order  to  obtain  the  damaged  elasto-plastic  stiffness  of  the  matrix  consfituent,  one 
needs  to  transform  eq.  (8)  from  the  undamaged  matrix  configuration  C®  to  the 
damaged  matrix  configuration  C®.  This  is  performed  through  the  material  time  differeo' 
tiation  of  eq.  (24)  together  with  strain  rate  counterpart  obtained  from  eqs  (49)  and 
(56),  such  that 

d€^  =  dJi/;:*7iS'  +  A/^*7dis.  (91) 

The  time  rate  of  the  matrix  damage  tensor  used  in  the  material  time  differentiation  of 
eq.  (24)  and  ite  inverse  used  in  eq.  (91)  may  be  expressed  as  shown  below  by  making 
use  of  eq.  (69) 


(83) 


(84) 


(85) 

(86) 

(87) 
#>e  easOy 

(88) 

# 


m  is  per* 
H}uations 
-followed 
com- 
ncorpor- 


dM 

(93) 

Making  use  of  eqs  (24)  and  (91)-(93),  one  obtains  the  resulting  elasto-plastic  stiffness  re¬ 
lation  for  the  damage  matrix  constituent 

rsfl' 

(94) 

where 

3A/™  8  ji/ 

Aim  _  V  m  cm  i  i  y  m  mnpq  m  r—m  cm 

—  0^"  ^  pqmrPmn  +  ^  ijkl  ~  '  0^_n)  ^  rtkl^pgcb^- 

(95) 

The  overall 
given  by 

response  of  the  composite  system  in  the  partial  effective  configuration. 

C,  is 

iSy  =  DykI  dl*/. 

(96) 

The  resulting  equation  for  D  is  given  by  [6] 

^{fkl  —  ^^pq^pqkl  ^ ^ypq^pqk!' 

(97) 

The  overall 
applying  the 

damage  response  of  the  composite  system  is  obtained  from  eq.  (96)  by 
interfacial  damage  effect  tensor  Using  the  following  relations 

45 y  —  5M  yi^u  4"  Jiu  do)u 

(98) 

and 

d(y  —  dM  +  M^i  dcu. 

(99) 

we  obtain  the  damage  elasto*plastic  constitutive  relation  including  both  the  local 
damages,  4)^  and  as  well  as  the  interfacial  damage  4>^.  Similarly,  the  rates  of  the 
debonding  damage  effect  tensor  used  in  eq.  (98)  and  its  inverse  used  in  eq.  (99)  are 
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given  as  follows  by  making  use  of  eq.  (87) 


,1  3A/1,  . 

(100) 

and 

g^  ^  pqfTS^rt- 

(101) 

Finally,  one  obtains  the  damage  elasto-plastic  constitutive  relation  including  both  the 
local  damages,  and  as  well  as  the  interfacial  damage  Making  use  of 

eqs  (96>-(101)  one  obtains 

d(Ty  =  Diju  dc*/. 

(102) 

where  the  damage  elasto-plastic  stiffness  of  the  material  is  given  by 

Dijkl  —  Of^ijbpqrsM 

(103) 

and 

^  _  1  »/d  A  imn  yi  F-l 

^Jkl  ~  ykl~  g^  rskl^  pqab^ab- 

(104) 

The  elastic  stiffness  for  the  damage  composite  £  is  given  such  that 

ffyt/  —  ^  pykJ' 

(105) 

where  the  elastic  stiffness  in  the  partial  effective  configuration  C  is  given 

as  follows: 

bijki  =  ^  Eyp^Ap^)./. 

(106) 

Thus, 

is 


The  c 
stre# 
equal: 


I 


when 

1 


# 


# 


COORDINATE  TRANSFORMATION 


The  three-dimensional  damage  elasto-plastic  constitutive  equation  for  single  lamina  refer¬ 
ring  to  the  principal  material  coordinate  system  has  been  introduced  in  eq.  (102).  The  general 
three-dimensional  constitutive  relation  of  a  composite  lamina  referring  to  the  off-axis  coordinate 
system  denoted  by  prime  can  be  obtained  from  eq.  (102)  by  coordinate  transformation. 
Here,  the  x-y  plane  coincides  with  the  xj-x:  plane,  and  the  angle  between  the  Xi  and  x  axis  is 
6.  The  stress  and  strain  vectors  in  those  two  coordinate  systems  are  related  by 


{do-}  =  IT]{d£r}' 

(de)»IT!{de}'. 

0 
0 
0 
0 

sin  0 
cos  6 

Substituting  eq.  (107)  into  eq.  (102),  we  obtain  the  relation 

{d<T}'  =  riT’[D]IT]{de}'. 


where  [T]  «  a  transformation  matrix  given  by 

r  cos  sin 

sin^S  cos^^ 

0  0 

cos  6  sin  0  —cos  0  sin  0 
0  0 


-2  cos  0  sin  0 
2  cos  0  sin  9 
0 

cos^0  -  sin^0 
0 


0 

0 

0 

0 

cos  9 
—sin  9 


(107) 


(108) 


(109) 


dinat 

coon 

sider 

ness 


Mak 

nate( 


I? 


(100) 


(101) 

;  both  the 
ng  use  of 

•  (102) 


(103) 


(104) 


(105) 


follows; 


(106) 


e  lamina  refer* 
2).  The  general 

coordinate 

vansformation. 

X,  and  X  axis  is 


(107) 


(108) 


J 
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Thus,  the  damage  clasto-plastic  stiflfness  matrix  referring  to  the  off-axis  coordinate  x-y-z  system 
is 


(110) 

The  constitutive  equation  for  the  plane  stress  problem  is  obtained  from  imposing  the  plane 
stress  conditions  ^  “  0  to  eq.  (109).  The  explicit  expression  of  constitutive 

equation  for  plane  stress  is  as  follows; 

dOjur 

'  - 

dOjcy 

where 

As 'x  2)3, 7^33' 

Z)*2  =  A2'-A3'xA27A3' 

A,  =  Ai '  ”  X  Ai  7A3' 

As  =  X  Dii'/Dii' 

Ai  =  Ai '  -  As '  X  2)31 7  As ' 

D22  =  2)22'  —  As'  X  Dii'/Dii 

As  =  A4 '  —  As '  X  2)34  7  As ' 

A2  *  A: '  —  As '  X  2)32  7A3 ' 

Daa' —  D43' X  DiA'/Dji'.  (112) 


'u  *'12  ^13 
^22  ^23 


ID 


tl 


2)?, 


A*,J 


d€*r 

6(yy 

2df 


xy 


(111) 


GROSS  DANUGE  ELASTO-PLASTIC  STIFFNESS 

The  elasto-plastic  damage  stiffness  tensor  for  a  single  lamina  in  its  principal  material  coor* 
dinate  system  has  been  presented  in  eq.  (103).  This  stiffness  tensor  is  transformed  to  the  loading 
coordinate  system  and  expressed  as  [D],c  in  matrix  form.  A  symmetric  stacking  of  plies  is  con* 
sidered  here  such  that  t  is  the  thickness  of  the  laminate  consisting  of  n  plies  and  is  the  thick* 
ness  of  the  kth  lamina.  The  average  stress  increment  is  expressed  as  follows  (in  vector  form); 

{d(r),v. »  [7£[D*l*^]{d<).  (113) 

Making  use  of  eq.  (113),  one  can  define  the  gross  damage  ela$to*plastic  stiffness  for  the  lami* 
Bated  composite  as  follows  in  matrix  form; 


(109) 


(114) 
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Making  use  of  the  assumption  of  constant  strain  through  the  laminate  thickness,  the  stresses  in 
each  lamina  are  calculated  as  follows: 


{d<T}*  =  [D  ]*{de). 


HNITE  ELEMENT  FORMULATION 

The  governing  equation  of  the  finite  element  method  can  be  derived  from  the  principle  of 
virtual  work  such  as 


I  ayS€ij  dK  =  I  q/Bui  dV  +  ^  tiBui  dA, 


where  Suj  is  a  field  of  virtual  displacements  that  is  compatible  with  applied  forces  and  Se^  is  the 
corresponding  field  of  compatible  virtual  strains  given  by 


,,  1  ra(5u,)  ,  a(a«,)l 


and  g,  and  t,  are  body  forces  and  surface  tractions,  respectively.  For  a  small  deformation  analy¬ 
sis,  we  have 

Ui^NijUj  (118) 

Bui^NyiBUj),  (119) 

where  Uj  is  the  displacement  of  nodal  points  and  Ny  is  the  displacement  interpolation  function 
or  the  shape  function. 

Substituting  eqs  (117)  and  (119)  into  eq.  (116),  one  obtains  the  equilibrium  equations  as  fol¬ 
lows: 


^  d  K  =  qiN>,  d  K  -I-  dA. 


One  finally  obtains  the  incremental  equilibrium  equations  by  differentiating  both  sides  of 
eq.  (120) 

[K]{dU)  =  {d?},  (121) 

where  {dU}  is  the  unknown  incremental  displacement  vector  of  the  nodal  points  and  {dP}  is  the 
corresponding  incremental  nodal  forces  given  by 


dP«  =  dqtNia  dF  -I- 1 dtiNu,  dA, 


where  dg,  is  the  incremental  body  force  and  dr*  is  the  incremental  surface  traction.  In  eq.  (121), 
[K]  is  the  stiffness  matrix  which  is  given  by 
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:  stresses  in 

(115) 

# 


The  incremental  equilibrium  eq.  (121)  expresses  the  equilibrium  between  Ae  internal  forces  {dF} 
(on  the  left-hand-side)  and  the  external  force  {dF}  (on  the  right-hand-side).  The  residual  force 
vector  {dF}  is  defined  by 

{dF}  =  {dF}-{dF).  (124) 


iprinciple  of 


•  (116) 


1 6tu  is  the 


In  a  damage  elastic— plastic  analysis,  because  of  the  nonlinear  relationship  between  Ae  stress 
and  strain,  the  equilibrium  eq.  (121)  is  a  nonlinear  equation  of  strains,  and  therefore,  is  a  non¬ 
linear  function  of  the  nodal  displacement.  Iterative  methods  are  usually  employed  to  solve 
eq.  (121)  for  displacements  conesponding  to  a  given  set  of  external  loads.  Moreover,  since  a 
damage  elasto-plastic  constitutive  relation  depends  on  deformation  history,  u  incremental 
analysis  following  an  actual  variation  of  external  forces  is  used  to  trace  the  variation  of  displace¬ 
ment,  strain,  stress  and  damage  along  with  the  external  forces. 

In  an  incremental  analysis,  the  total  load  {F}  acting  on  a  structure  is  added  in  increments 
step  by  step.  At  the  (n  +  l)th  step,  the  load  can  be  expressed  as 

"+•{/*)  ="(J’)+'^MdF},  (125) 

where  the  left  superscript  n  indicates  the  nth  incremental  step.  Assuming  that  the  solution  at  the 
nth  step,  "{u},  "{ff},  "{c}  and  "{d>}  is  known,  and  at  the  (n  +  l)th  step,  one  obuins  the  follow- 


(117) 


ing,  corresponding  to  the  load  increment  {dF}, 


# 

"+’{«)  ="  {«}  +  ld«) 

(126) 

ion  analy- 

"+'{a)  ="  {a}-l-{da} 

(127) 

(118) 

"^'{€)="{f}  +  [d€} 

(128) 

(119) 

+  r  =  m,  f,  d. 

(129) 

1  function 

^ns  as  fol- 

SOLUTION 

(120) 

^sides  of 
(121) 
is  the 

(122) 
\  (121), 


A  full  Newton-Raphson  method  is  used  in  this  work  to  solve  the  system  of  nonlinear 
equations  that  arise  from  the  equilibrium  equations.  A  brief  description  of  the  method  is  given  by 
Voyiadjis  [13].  The  incremental  analysis  technique  described  in  this  chapter  is  successfully  im¬ 
plemented  into  the  finite  element  program  NDA  (Nonlinear  Damage  Analysis)  using  the  above 
described  iterative  method.  The  steps  involved  in  the  process  of  solving  are  briefly  described  below. 

(1)  INCREMENT:  Loop  for  each  load  increment 

(1)  Calculate  the  load  or  applied  displacement  increment  for  the  current  incremental  step 
or  input  the  load/applied  displacement  increment. 

(2)  ITERATE:  Loop  for  full  Newton-Raphson  iteration: 

(1)  Compute  the  residual  load  vector  for  this  iteration  subtracting  the  equilibrium 
load  from  the  load  computed  for  the  increment. 

(2)  Rotate  the  appropriate  loads  and  applied  displacements  such  that  the  degrees 
of  freedom  at  the  skew  boundary  (a  boundary  condition  that  is  not  along  the 
global  coordinate  system)  are  normal  and  tangential  to  the  skew  boundary. 

(3)  Assemble  the  stiffness  matrices  and  find  the  equivalent  loads  for  the  applied 
incremental  displacements.  Since  explicit  integration  is  difficult,  Gaussian 
points  are  used  to  evaluate  the  above  integrals. 

(4)  Solve  for  the  incremental  displacements  using  a  linear  solver. 
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(5)  Add  the  solved  iterative  incremental  displacements  to  the  applied  incremental 
displacements  to  obtain  the  complete  iterative  incremental  displacements. 

(6)  Rotate  back  the  complete  iterative  incremental  displacements  at  the  skew 
boundaries  to  the  global  coordinate  system. 

(7)  Cumulate  the  complete  iterative  incremental  displacements  to  the  total  incre- 
menul  displacements. 

(8)  Find  the  stresses  due  to  the  iterative  incremental  displacements.  From  the 
iterative  deformation  gradient  and  the  stresses  updated,  compute  the  updated 
constitutive  matrix  D.  From  the  total  incremental  displacements  accumulated 
so  far  and  the  D  matrix,  calculate  the  equilibrium  load  vector. 

(9)  Check  if  the  convergence  of  solution  is  met  using  a  particular  convergence  cri¬ 
terion.  If  convergence  has  not  occurred,  go  back  to  the  step  ITERATE. 

(10)  If  divergence  occurs  according  to  the  convergence  criterion,  then  reduce  the 
load  increment  appropriately  as  specified  by  the  user  and  start  the  iterative 
solution  over  again  for  that  load  increment. 

(1 1)  If  divergence  occurs  for  a  load  increment  that  has  been  reduced  “m”  times 
(specified  by  the  user),  then  report  “convergence  not  met”  and  leave  the  sol¬ 
ution  phase. 

(12)  If  convergence  has  occurred,  then  perform  the  following  operations  before 
going  for  the  next  increment. 

(1)  Update  the  nodal  positions  by  adding  the  currently  obtained  incre¬ 
mental  displacements. 

(2)  Transform  the  quantities  pertaining  to  the  material  property  to  the 
present  configuration. 

(3)  Print  out  the  appropriate  quantities  pertaining  to  the  converged  incre¬ 
ment  according  to  the  user’s  specifications. 

(4)  If  the  total  load  is  not  reached,  go  back  to  the  step  INCREMENT. 


# 


(1 
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Fig.  2.  Finite  element  meshes. 


STRESS  AND  DAMAGE  COMPUTATIONS 


(1)  Step  1.  Retrieve  ay,  a'y,  4>y.  Retrieve  also  the  information  whether  the  previous  loading 
was  a  damage  loading  or  not  (IDAMG)  and  plastic  loading  or  not  (lYILD). 

(1)  If  IDAMG  =  0  when  retrieved,  then  evaluate  the  incremental  elastic-predictor 
stress  ffj  assuming  that  the  loading  is  elastic.  Use  the  undamaged  elastic  stiffness 
matrix  for  the  calculation  (da^j^Eykt  dtiu). 

(2)  If  IDAMG  #  0  when  retrieved,  use  (dag  *  Eyki  dc*/). 

(3)  Calculate  the  incremental  elastic-predictor  stress  of  matrix  constituent 

(4)  Check  if  the  predicted  stress  state  of  matrix  constituent  is  inside  the  yield  surface 
or  not. 


Table  1.  Material  properties 


Matrix  {Ti-i4Al-21Nb) 

Fiber  (SiQ 

Modulus 

8x10*  MPa 

41x10*  MPa 

Poisson*!  ratio 

0.30 

0.22 

Initial  volume  fraction 

0.65 

0.35 

Yielding  stress  ^ 

360  MPa 

Kinematic  hardening  parameter  b 

90  MPa 

Table  2.  Local  damage  parameters 


Matrix  damage 

Fiber  damage 

Interfacial  damage 

0.08 

0.06 

0.075 

'll 

0.08 

0.06 

0.065 

•la 

0.08 

0.06 

0.065 

'll 

0.65 

0.55 

0.55 

Si 

0.65 

0.55 

0.70 

0.65 

0.55 

0.70 

Si 

m 

0.003 
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Fig.  3.  (a)  Stress-strain  curves  of  I±45]s  layup,  (b)  Stress-strain  curves  of  t0/90)s  layup. 

(5)  If  the  stress  state  of  matrix  constituent  is  inside  the  yield  surface  then: 

(1)  Assign  elastic  stiffness  to  the  constitutive  stiffness  and  the  predictor  stress 
increment  to  the  actual  computed  stress  increment. 

(2)  Set  lYILD  =  0  indicating  the  elastic  loading  has  uken  place. 

(3)  Exit  to  Step  2.  Otherwise,  go  to  the  next  step. 

(4)  Set  lYILD  =  I.  then: 

(1)  Calculate  the  elasto-plastic  stiffness  D  (when  IDAMG  =  1)  or  .D 
(when  IDAMG  =  0). 

(2)  Update  the  quantities  Oy,  Oy,  a™. 


08  10  1S  20  2S303S  0S10  1S  20  2SM3S  0S10  1S  20  2S  30  35 

«a  [*<51,  {145],  1145], 


Fig.  4.  Strain  contoun  for  [±45]$  layup  (in  •/•). 


Fig.  5.  Stnio  contours  for  (0/90]s  byup  (in  */•). 


(2)  Step  2 

(1)  Check  the  damage  criteria  using  the  updated  quantity  <7^-. 

(2)  If  damage  criteria  g'  <  0,  then  IDAMG  “  0.  Exit  from  the  routines. 

(3)  If  damage  criteria  g^>0,  then  IDAMG  *  1.  Calculate  the  damage 
increment  d(^'  and  update  damage  quantity 

(4)  Store  the  updated  quantities  in  a  file. 


APPLICATION  TO  THE  DOG-BONE  SHAPED  SPECIMEN  AND  THE  CENTER- 
CRACKED  LAMINATED  PLATES 

The  finite  element  method  is  used  for  solving  a  dog-bone  shaped  specimen  and  a  center- 
cracked  laminate  plate  shown  in  Fig.  1  that  is  subjected  to  inplane  tension.  Due  to  symmetry  in 
geometry  and  loading  as  shown  in  Fig.  1,  one-quarter  of  the  plate  needs  to  be  analyzed.  Two- 
dimensional  plane  stress  analysis  rather  than  three-dimensional  analysis  is  used  here  since  the 
thickness  of  plate  is  much  smaller  than  the  other  dimensions.  Applying  the  appropriate  bound¬ 
ary  conditions  for  the  symmetry,  both  one-quarter  of  the  center-cracked  laminate  plate  and  the 
dog-bone  shaped  specimen  are  discretized  using  plane  stress  finite  elements.  The  finite  element 
meshes  chosen  for  analyzing  the  problems  are  shown  in  Fig.  2.  The  fom-noded  quadrilateral  el¬ 
ement  is  used  in  both  finite  element  analyses. 

Two  types  of  laminate  layups  (±45),  and  (0/90),  each  consisting  of  four  plies  are  used  here. 
The  thickness  of  each  ply  is  equal  to  0.254  mm.  Since  both  layups  are  symmetric,  no  curvature 
is  assumed.  Hence,  the  strain  through  the  plate  thickness  is  assumed  to  be  the  same.  The  ma- 


Fig.  6.  Stress  eontoun  for  I±45]s  Uyup  (units  are  in  MPa). 


terial  properties  and  damage  parameters  using  the  proposed  constitutive  model  are  listed  in  ij  Coi 

Tables  1  and  2,  respectively.  h 

The  following  convergence  criterion  is  used  in  this  analysis  which  is  based  on  the  incremen-  ’’ 

tal  internal  energy  for  each  iteration  in  that  incremental  loadingtl4].  It  represents  the  amount  of  ' 

work  done  by  the  out-of-balance  loads  on  the  displacement  increments.  Comparison  is  made  th? 

with  the  initial  internal  energy  increment  to  determine  whether  or  not  convergence  has  occurred.  ten 
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Fig.  7.  Stress  contours  for  [0/90]$  Uyup  (uniu  are  in  MPa). 

Convergence  is  assumed  to  occur  if  for  an  energy  tolerance  £e,  the  following  condition  is  met: 

/>«-!))  <  -"/)),  (130) 

where  is  the  incremental  displacement  residual  at  the  (Oth  iteration,  T  * 

the  out-of-balance  force  vector  at  (f  - 1)  iteration  and  (At^V  */?-"/))  «  the  internal  energy 

term  for  the  (Oth  iteration  in  the  (n  +  l)th  increment.  Divergence  is  assumed  to  occur  if  the 
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Fig.  8.  Comparison  of  damage  elasto-plastic  analysis  with  elastic  analysis  of  stress  Oyy  contours  around 
the  crack  tip  for  (±45)$  layup  (units  are  in  MPa). 
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out-of-balancc  internal  energy  for  the  (/  -  l)th  iteration  is  greater  than  the  out*of-balance  in¬ 
ternal  energy  for  the  (0th  iteration. 

The  load  is  incremented  with  uniform  load  increments  of  5  MPa  until  the  principal  maxi¬ 
mum  local  damage  value  4>p  reaches  1 .0  ((pp  ^  1 .0).  The  principal  maximum  local  damage  value 
4>p  is  given  by: 


r  =  m,  f ,  d. 


(131) 


Consequently',  material  failure  at  integration  point  is  assumed  when  k  1.  The  principal 
damage  value  of  the  integration  point  in  all  elements  is  monitored  at  each  load  increment  since 
it  is  used  to  determine  the  onset  of  macro<rack  initiation  of  the  material. 

The  dog-bone  shaped  specimen  failed  when  the  final  load  of  270  MPa  was  reached  for  the 
(±45)s  layup  and  480  MPa  for  the  (0/90)$  layup.  These  failure  loads  are  close  to  the  experimen¬ 
tal  failure  loads  276  MPa  for  the  (±45)s  layup  and  483  MPa  for  the  (0/90)s  layup  [15].  The  ma¬ 
terial  failure  for  the  center-cracked  specimen  occurs  at  the  front  of  the  crack  tip  when  the  final 
load  of  80  MPa  is  reached  for  the  (±45)s  layup  plate  and  120  MPa  for  the  (0/90)s  layup  plate. 


DISCUSSION  OF  THE  RESULTS 

The  stress-strain  curves  from  both  the  finite  element  analyses  and  experiments  of  the  two 
types  of  layups  of  the  dog-bone  shaped  specimens  are  shown  in  Fig.  3.  Good  correlation  is 
shown  between  the  finite  element  analysis  results  and  the  experimental  data  obtained  by 
Voyiadjis  and  Vensontl5]. 
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Fig.  9.  Comparison  of  damage  e)asto*plastic  analysis  with  elastic  analysis  of  stress  Cyy  contours  around 
the  crack  tip  for  (0/90j^  layup  (units  are  in  MPa). 

Strain  contours  for  the  (±45)s  layup  and  (0/90)s  layup  of  the  center  cracked  plates  are 
shown  in  Figs  4  and  5,  respectively.  Since  the  two  types  of  layups  are  symmetric,  the  strains  in 
each  laminae  of  the  layup  are  the  same.  However,  the  stress  and  damage  distributions  are  differ¬ 
ent  for  each  laminae  of  the  layup  since  each  laminae  has  a  different  stiffness.  Stress  contours  for 
each  laminae  are  indicated  in  Fig.  6  for  the  (±45)s  layup  and  Fig.  7  for  the  (0/90)$  layup.  In 
Figs  8  and  9,  comparison  is  made  between  the  damage  analysis  and  the  elastic  analysis  for  the 
stress  Cyy  contours  around  the  crack  tip.  The  damage  analysis  shows  considerable  stress  re¬ 
duction  due  to  the  damage  around  the  crack  tip.  The  stress  Oyy  at  the  front  of  the  crack  tip  as 
obtained  from  the  elastic  solution  is  higher  than  that  of  the  material  strength  of  the  layup. 
However,  in  the  damage  elasto-plastic  analysis,  the  stresses  are  reduced  such  that  they  are  close 
to  those  of  the  material  strength.  The  Pyy  stress  reductions  at  the  front  of  the  crack  tip  are  more 
than  50%  for  [±45],  40%  for  [0]  ply  and  80%  for  [90]  ply.  Stress  redistributions  are  clearly  indi¬ 
cated  in  Figs  8  and  9.  Primarily  due  to  the  stress  reduction  around  the  crack  tip,  the  stress  is 
therefore  transferred  to  the  outer  portion  away  from  the  crack  tip.  This  is  clearly  indicated  in 
Fig.  9  where  the  stress  reduction  at  the  90^  ply  is  primarily  due  to  considerable  interfacial 
damage. 

The  local  damage  contours  around  the  crack  tip  are  shown  in  Figs  10-13,  for  the  failure 
loads  in  the  case  of  [+45],  [—45],  [0],  and  [90]  ply,  respectively.  For  the  [±45]  layups,  all  types  of 
damage  such  as  matrix,  fiber  and  interfacial  are  developed.  Fiber  damage  is  considerably  more 
spread  in  the  [0]  ply  than  the  interfacial  damage.  On  the  other  hand,  interfacial  damage  is  more 
pronounced  with  matrix  damage  for  the  [90]  ply.  However,  fiber  damage  is  much  less  developed 
in  the  case  of  the  [90]  ply.  This  is  in  line  with  the  experimental  results  obtained  by  Voyiadjis 
and  Venson[15]. 
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Fig.  10.  Damage  contours  around  crack  tip  at  the  failure  load  for  [■*■45]  lamine 


SUMMARY  AND  CONCLUSIONS 

TTie  proposed  constitutive  model  is  implemented  numerically  using  the  finite  element 
method.  The  model  is  used  to  analyze  the  dog-bone  shaped  specimens  and  the  center-cracked 
laminated  plates  subjected  to  inplane  tensile  forces.  Very  good  correlations  are  demonstrated 
between  the  numerical  results  obtained  using  the  proposed  theories  and  the  experimental  results 
for  uniaxial  tension.  The  stress  and  damage  contours  in  the  case  of  the  center  cracked  plate 
show  that  stress  redistributions  and  damage  are  qualitatively  in  line  with  the  physics  of  defor¬ 
mation.  The  analysis  presented  here  allows  the  separate  quantification  of  the  different  types  of 
damages  such  as  matrix,  fiber  or  debonding. 
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Fig.  11.  Damage  contoun  around  crack  tip  at  the  failure  load  for  [-45}  lamina. 
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Fig.  12.  Damage  contoun  around  crack  tip  at  the  failure  load  for  (0]  lamina. 
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Fig.  13.  Damage  contoun  around  crack  tip  at  the  failure  load  for  [90]  lamina. 

The  authors  are  currently  working  on  damage  due  to  delamination  which  will  be  introduced 
into  the  proposed  model  in  future  work.  In  order  to  capture  delamination  due  to  interlamina 
stresses,  a  three^Jimensional,  finite  element  analysis  will  be  performed. 
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APPENDIX  A 

Mori-Tanaka's  strain  and  stress  concentration  tensors 

The  expressions  for  the  elastic  stress  and  strain  concentration  facton  given  here  are  based  on  the  Mori-Tanaka  method. 
In  the  recent  paper  by  Chen  et  al.[\€\,  the  expressions  for  the  elastic  strain  concentration  factors  A'  and  the  clastic  stress 
concentration  factors  S'  are  given  by 


A^ 

r»  m,  f 

(AI) 

r  =  m,  f. 

(A2) 

(A3) 

(A4) 

The  tensors  and  are  termed  the  partial  concentration  facton  for  strain  and  stress,  and  arc  expressed  in  the  follow¬ 
ing  form: 

(A5) 

(A6) 

(a?) 


where  «nd  £“  are  the  elMtic  stiffness  tensors  of  the  fiber  and  matrix,  respectively.  The  tensors  of  ^  and  C  depend 
only  on  the  shape  of  the  inclusion  and  the  elastic  moduli  of  the  surrounding  matrix.  For  example,  for  an  inclusion  in  the 
shape  of  a  circular  cybnder  in  isouopic  matrix,  the  tensor  P  written  in  matrix  form  (6x6  array)  is  given  by 
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where  £®  is  Young’s  modulus  of  the  matrix,  i®  is  the  Poisson  ratio  of  the  matrix  and  is  the  shear  modulus  of  the 
matrix.  The  tensor  C  in  cq.  (16)  is  given  by 
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Based  on  the  concept  of  the  effective  stress  and  on  the  description  of  anisotropic 
damage  deformation  within  the  framework  of  continuum  damage  mechanics,  a  fourth 
order  damage  effective  tensor  is  properly  defined.  For  a  general  state  of  deformation 
and  damage,  it  is  seen  that  the  effective  stress  tensor  is  usually  asymmetric.  Its 
symmetrization  is  necessary  for  a  continuum  theory  to  be  valid  in  the  classical  sense. 
In  order  to  transform  the  current  stress  tensor  to  a  symmetric  effective  stress  tensor, 
a  fourth  order  damage  effect  tensor  should  be  defined  such  that  it  follows  the  rules 
of  tensor  algebra  and  maintains  a  physical  description  of  damage.  Moreover,  an 
explicit  expression  of  the  damage  effect  tensor  is  of  particular  importance  in  order 
to  obtain  the  constitutive  relation  in  the  damaged  material.  The  damage  effect  tensor 
in  this  work  is  explicitly  characterized  in  terms  of  a  kinematic  measure  of  damage 
through  a  second-order  damage  tensor.  In  this  work,  tensorial  forms  are  used  for 
the  derivation  of  such  a  linear  transformation  tensor  which  is  then  converted  to  a 
matrix  form. 


Introduction 

In  1958.  Kachanov  (Kachanov,  1958)  introduced  the  concept 
of  effective  stress  in  damaged  materials.  This  pioneering  work 
started  the  subject  that  is  now  known  as  continuum  damage 
mechanics.  Research  in  this  area  has  steadily  grow'n  and  reached 
a  stage  that  warrants  its  use  in  today’s  engineering  applications. 
Continuum  damage  mechanics  is  now  widely  used  in  different 
areas  including  brittle  (Krajcinovic,  1983;  Krajcinovic  and 
Foneska,  1981)  and  ductile  failure  (Lemaitre,  1985,  1986).  In 
this  theory,  a  continuous  damage  variable  is  defined  and  used 
to  represent  degradation  of  the  material  which  reflects  various 
types  of  damage  at  the  microscale  level  like  nucleation  and 
growth  of  voids,  cavities,  microcrack,  and  other  microscopic 
defects. 

In  continuum  damage  mechanics,  the  effective  stress  tensor 
is  usually  not  symmetric.  This  leads  to  a  complicated  theory  of 
damage  mechanics  involving  micropolar  media  and  the  Cosserat 
continuum.  Therefore,  to  avoid  such  a  theory,  symmetrization 
of  the  effective  stress  tensor  is  used  to  formulate  a  continuum 
damage  theory  in  the  classical  sense.  Several  methods  used  in 
order  to  symmetrize  the  effective  stress  tensor  are  proposed  by 
Lee  et  al.  ( 1986)  and  Sidoroff  ( 1979).  A  linear  transformation 
tensor,  defined  as  a  fourth-order  damage  effect  tensor,  is  pro¬ 
posed  by  Sidoroff  (1979);  however,  no  explicit  form  of  this 
tensor  is  given.  Moreover,  damage  tensor  of  higher  order  than 
two  may  fail  to  convey  the  physical  meaning  of  damage.  In 
addition,  the  works  of  Lee  et  al.  ( 1986)  and  Sidoroff  ( 1979)  are 
confined  to  two-dimensional  problems.  Furthermore,  no  explicit 
expressions  are  derived  for  ^e  fourth-order  linear  transforma¬ 
tion  tensors  for  the  general  anisotropic  damage  behavior  of 
three-dimensional  problems. 

In  this  work,  continuum  damage  mechanics  will  be  reviewed 
based  on  the  concept  of  effective  stress.  The  effective  stress 
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is  defined  as  the  stress  acting  on  a  hypothetical  undamaged 
configuration  that  produces  the  same  clastic  strain  or  elastic 
strain  energy  as  the  actual  state  of  stress  acting  on  the  current 
damaged  configuration  based  on  the  equivalence  hypothesis. 
This  equivalence  statement  is  known  as  the  hypothesis  of  elastic 
strain  equivalence  or  the  hypothesis  of  elastic  energy  equiva¬ 
lence.  Using  the  definition  of  the  effective  stress  and  the  hypoth¬ 
esis  of  elastic  energy  equivalence,  one  can  solve  for  variable  in 
the  hypothetical  undamaged  configuration,  such  as  the  effective 
strain.  However,  strain  in  the  hypothetical  undamaged  configu¬ 
ration  is  equal  to  that  in  the  current  damaged  configuration 
under  the  hypothesis  of  elastic  strain  equivalence.  For  a  detailed 
review  of  the  principles  of  continuum  damage  mechanics  as 
used  in  this  work,  the  reader  is  refened  to  the  works  of  Kacha¬ 
nov  ( 1958),  Lemaitre  ( 1985,  1986),  Krajcinovic  ( 1985),  Cha- 
boche  (1981,  1988a,  b),  Murakami  (1988),  Sidoroff  (1979, 
1980),  and  Voyiadjis  and  Kattan  ( 1992). 

In  a  general  state  of  deformation  and  damage,  the  effective 
stress  tensor  et  is  related  to  the  stress  tensor  a  by  the  following 
linear  transformation: 

^ij  ~  MikjiCTu  (1) 

where  <r  is  the  Cauchy  stress  tensor  and  M  is  a  fourth-order 
linear  transformation  operator  called  the  damage  effect  tensor. 
Depending  on  the  form  used  for  M,  it  is  very  clear  from  Eq. 

( 1 )  that  the  effective  stress  tensor  et  is  generally  not  symmetric. 
Using  a  nonsymmetric  effective  stress  tensor  as  given  by  Eq. 

( 1 )  to  formulate  a  constitutive  model  will  result  in  the  introduc¬ 
tion  of  the  Cosserat  and  micropolar  continua.  However,  the  use 
of  such  complicated  mechanics  can  be  easily  avoided  if  the 
proper  fourth-order  linear  transformation  tensor  is  formulated 
in  order  to  symmetrize  the  effective  stress  tensor.  Such  a  linear 
transformation  tensor  called  the  damage  effect  tensor  is  obtained 
in  the  literature  (Lee  et  al.,  1986;  Sidoroff,  1979)  using  symme¬ 
trization  methods.  However,  it  lacks  a  systematic  and  consistent 
approach.  It  is  the  aim  of  this  work  to  provide  a  solid  basis  for 
such  transformation  of  the  second-order  stress  tensor  and  its 
justification  for  the  symmetrization.  Three  different  formula¬ 
tions  for  symmetrization  of  the  effective  stress  tensors  proposed 
by  Lee  et  ( 1986)  and  Sidoroff  ( 1979)  are  described  below. 
The  effective  stress  tensor  is  symmetrized  using  the  following 
laws: 
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^ij  ~  (^lA  ”  <t>tk)  “  <#>;/)  (2) 

Vij  =  \[crA6,,  -  +  (6, A  -  (3) 

~  —  ^A;)  +  (^;A  ”  <t>jk)^ki]^  (4) 

Other  forms  for  symmeirization  are  defined  by  many  researchers 
through  the  fourth-order  damage  effect  tensor  such  as  that  by 
Cordebois  and  Sidoroff  (1982) 

=  (^.A  -  (5) 

by  Murakami  and  Ohno  (1980) 

^tkjl  ”  ^  *“  ^kj)~^ 

+  (<5y  -  (i>>jy'^bu  +  (6,/  -  <i>,iy^bkj],  (6) 

by  Betten  (1983) 

^ikji  =  ~  <^A/)  * 

^(S„-4>,y^(6,j-<}>,X')]^  (7) 

and  by  Lu  and  Chow  (1990) 

Miiji  =  5[cxp((f>,/2)  exp(<^ij/2) 

+  exp(<f>,,/2)  exp(<;>i/2)].  (8) 

in  the  above  equations  is  a  damage  tensor  characterized  by 
a  second-order  symmetric  tensor  and  is  given  by  (Murakami, 
1983) 

3 

<i>ij  =  X  (no  sum  in  it)  (9) 

A=1 

where  n*  is  an  eigenvector  corresponding  to  the  eigenvalue. 
<#)a,  of  the  damage  tensor.  Voyiadjis  and  Venson  (1995) 
quantified  the  physical  values  of  the  eigenvalues  (A-  =  L  2, 
3 )  and  the  second-order  damage  tensor  4>  for  the  unidirectional 
fibrous  composite  by  measuring  the  crack  density  with  the  as¬ 
sumption  that  one  of  the  eigendirections  of  damage  tensor  coin¬ 
cides  with  the  fiber  direction.  This  introduces  a  distinct  kine¬ 
matic  measure  of  damage  which  is  complimentaiy'  to  the  defor¬ 
mation  kinematic  measure  of  strain.  A  thermodynamically 
consistent  evolution  equation  for  damage  tensor  4>  together  with 
a  generalized  thermodynamic  force  conjugate,  Y,  to  the  damage 
tensor  is  presented  in  the  paper  by  Voyiadjis  and  Park  ( 1995). 

Numerous  fourth-order  damage  effect  tensors  are  defined  by 
using  the  symmetrizatipn  laws  indicated  above.  However,  only 
the  one  by  Cordebois  and  Sidoroff  (1982)  may  be  obtained 
from  the  symmeirization  procedure  given  by  Eq.  (2).  One  can¬ 
not  deduce  explicitly  the  fourth-order  damage  effect  tensor  M 
from  the  remaining  proposed  procedures.  In  the  case  of  Corde¬ 
bois  and  Sidoroff  (1982)  it  is  impossible  to  get  the  explicit 
form  of  the  square  root  of  the  second-order  tensor  in  Eq.  (5). 
Alternatively,  the  damage  effect  tensor  using  the  fourth-order 
damage  tensor  is  defined  by  Chaboche  (1979)  as  follow's: 

M,kjl  =  (likjl  “  *Pikjl)  ^  (10) 

where  /,aj/  is  a  fourth-order  identity  tensor  and  is  given  by 

Iikji  =  +  ^ii^kj)'  (11) 

However,  it  is  not  easy  to  characterize  physically  the  fourth- 
order  damage  tensor  ipijki  rather  than  the  second-order  damage 
tensor  4>,j.  For  the  case  of  isotropic  damage,  the  fourth-order 
damage  tensor  is  defined  by  Ju  ( 1990)  as  follows: 

^ikji  =  d]6ik6ji  -f  djlikjf  (12) 

where  and  d2  are  scalar  (dependent  or  indep)endent )  damage 
variables.  Using  the  second-order  anisotropic  damage  tensor  4>,j 


in  the  damage  effect  tensors  given  by  Eqs.  (5),  (6),  (7),  and 
(8)  one  may  lose  the  physical  view  of  the  net  stress  tensor  due 
to  the  presence  of  the  off  diagonal  elements  of  the  damage 
tensor  4>,j.  In  order  to  avoid  this  problem,  the  principal  damage 
tensor  rather  than  the  second-order  damage  tensor  is  used  in 
conjunction  the  damage  effect  tensor.  However,  the  eigendirec¬ 
tions  of  the  damage  tensor  do  not  coincide  with  the  stress  tensor 
or  eigendirection  of  the  stress  tensor.  Since  the  damage  tensor 
^  always  has  three  orthogonal  principal  directions  n^(A  =  1, 
2,  3)  and  the  three  corresponding  principal  values  <#>a(A  =  1,2, 
3),  Eqs.  (2),  (3),  and  (4)  can  be  expressed  as  follows  in  the 
coordinates  that  coincide  with  the  three  orthogonal  principal 
directions  of  the  damage  tensor: 

Pm.  =  (6mr  -  ^mrr''^Op,(6r^  “  (13) 

Pm.  =  +  ( <5.r  “  (14) 

&m.  =  ^lPmr(6r.  "  ^r.)  +  (6.r  "  ^.r)Prm]  (15) 
where  is  a  principle  damage  tensor  that  is  given  by 

^ij  =  hikbji(l>u 

0  ol 

=  0  ^2  0  (16) 

.  0  0  <^>3j 

and  the  second-order  transformation  tensor  b  is  given  by 


n2 

2 

rr] 

nl 

riy 

m3 

m3 

^2 

m3 

This  transformation  tensor  called  the  proper  orthogonal  tensor 
requires  that 

(18) 

The  effective  stress  tensor  in  the  principal  damage  direction 
coordinates  system  is  given  by 

9m.  =  (19) 

Similarly,  the  stress  tensor  in  the  principal  damage  direction 
coordinates  system  is  given  by 

=  bpkb^iCTu.  (20) 

Using  the  principle  damage  direction  coordinate  system,  Eq. 
( 1 )  is  given  by 

^ mr>  ^dffipnqO pq  (21  ) 

The  fourth-order  damage  effect  tensors  given  by  Eqs.  (5),  (6), 
(7),  and  (8)  should  be  now  expressed  as  follows: 

=  (^m.  -  (22) 

^mpfuj  “  l\^mn{bpq  —  ^pq)  ~  ^pn) 

+  (<5^  --  ^mI^y'bpq  +  (<5^  -  ^mqy'bp„]  (23) 

^mpruj  ^m.)  “  ^pq) 

^  4>^qyH6p„-  ^pny^)]  (24) 

^mpnq  ““  J  )  ®^P(^p9^2) 

+  exp((^^/2)cxp((#i;,„/2)],  (25) 

respectively.  These  tensors  are  termed  the  principal  damage 
effect  tensors. 
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Fourth-Order  Anisotropic  Damage  Effect  Tensor 

The  explicit  representation  of  the  fourth-order  damage  effect 
tensor  M  using  the  second-order  damage  tensor  is  of  particu¬ 
lar  importance  in  the  constitutive  modeling  of  damage  mechan¬ 
ics.  However,  it  is  impossible  to  use  the  damage  tensor  rather 
than  the  principal  damage  tensor  ^  directly  in  the  formulation. 
Therefore  the  damage  effect  tensor  M  in  Eq,  (1)  should  be 
obtained  from  Eq.  (21 )  using  coordinate  transformation. 

Substituting  Eqs.  (19)  and  (20)  into  Eq.  (21),  one  obtains 
the  following  relation: 

^ij  ^m^njPpJ^ql^mpnq^ kl-  (26) 

Therefore  the  fourth-order  tensor  M  in  Eq.  (1)  is  given  as 
follows: 


In  order  to  obtain  the  matrix  form  of  the  damage  eucci 
the  stress  tensors  dy  and  7*;  in  Eq.  ( 1 )  are  converted  into  vector 
form  as  follows: 

{ a)  =  { dn  0*22  c^33  <3^12  cr23  (Ti3  }  ^  (33) 

(7)  =  {2Tii  ^22  ^33  2^12  ^23  (34) 

Using  the  notation  of  Eqs,  (33)  and  (34),  Eq.  (1)  is  now 
represented  in  matrix  notation  as  follows: 

{^}  =  (35) 

where  [A/]  is  the  6  x  6  matrix  representation  of  the  fourth- 
order  tensor  M.  The  explicit  form  of  the  matrix  [M]  is  given 
as  follows: 


2di2d]3 

2^22^23 

2^32033 

(dl2d23  +  ^13^22) 
(^22^33  +  ^23^32) 
(^i2^33  +  fll3<33l) 


ai2<3i2 

fl)3ai3 

20]\0i2 

^2)^21 

022021 

^23^23 

2O21O22 

^^31^31 

^32^32 

<233^^33 

2O31O32 

^11^^21 

012022 

^23^23 

(0i]022  +  fll2^2l) 

^21^31 

022^32 

^223^33 

(O21O32  +  fl23^3l) 

^11^31 

O12O32 

13^33 

(011^32  +  Oj303j) 

2dndl3 

2^21^23 

2d3jd33 

(^11^23  +  fll3^2l) 
(021^33  4-  023^31) 
(^11^33  +  fll3^3l) 


(36) 


A/,*y/  =  b„ibnjbpkb^iM„ 


(27) 


It  is  clear  that  the  fourth-order  damage  effect  tensors  presented 
by  Eqs.  (5),  (6),  (7),  and  (8)  differ  from  the  damage  effect 
tensor  obtained  by  Eq.  (27).  Therefore  the  fourth-order  damage 
effect  tensor  presented  by  the  Eqs.  ( 5 ) ,  ( 6 ) ,  ( 7 ) ,  and  ( 8 )  should 
be  expressed  in  the  princpal  damage  direction  coordinate  system 
using  the  principal  damage  tensor  ^ . 

One  of  the  explicit  expression  for  the  fourth-order  damage 
effect  tensors  using  the  principal  damage  effect  tensor  given  by 
Eq.  (22)  is  presented  here.  The  principal  damage  effect  tensor 
given  by  Eq.  (22)  can  be  written  as  follows: 

=  (28) 


where  the  second-order  tensor  a  is  given  by 
1 


\1  ~ 
0 

0 


V  1  ~  <#>2 


1 


VI  -</>2  J 


(29) 


Substituting  Eq.  (28)  into  Eq.  (27),  one  obtains  the  following 
relation: 


Mikjl  =  b^,bnjbpJ>q4mf4nq 


=  (30) 

Using  Eq.  (30),  a  second-order  tensor  a  is  defined  as  follows: 

Oik  -  b^ibpi4mp-  (31) 

The  matrix  form  of  Eq,  (31 )  is  as  follows: 


It  is  clear  that  the  various  explicit  representation  forms  of  the 
fourth-order  damage  effect  tensors  given  by  Eqs.  (5), (6), (7), 
and  ( 8 )  using  the  second  damage  tensor  are  violated  by  the 
tensor  transformation  law. 

Matrix  Forms  of  the  Damage  Effect  Tensors  for  Two- 
Dimensional  Problems 

The  explicit  matrix  forms  of  the  damage  effect  tensors  for 
two-dimensional  problems  are  presented  in  this  section.  For  the 
shake  of  simplicity,  a  plane  state  of  damage,  =  0  or  <^>33  = 
(^>23  =  </>32  =  <^>13  =  <#>3]  =  0,  is  assumed  for  both  plane  stress 
and  plane  strain  problems.  Similar  stress  vector  of  Eqs.  (33) 
and  (34),  for  the  two-dimensional  problems  are  given  by 

(a)  =  {ff,,  (7:2  (712)'’  (37) 

{?)  =  { (38) 

Using  the  notation  of  Eqs.  (37)  and  (38),  one  obtains  Eq.  (21 ) 
represented  in  matrix  notation  as  follows: 

{^)  =  [M][d]  (39) 

where  the  3  X  3  matrix,  [M],  termed  principal  damage  effect 
matrix  is  given  by 


Afnn 

^M\l 

Afin2 

A/i2I1 

A/2121 

A/2222 

^2122 

+ 

A/222I 

(40) 

^1121 

A/1222 

M\\22 

+ 

A/j221_ 

The  principal  damage  effect  matrices  corresponding  to  Eqs. 
(22),  (23),  (24),  and  (25)  are  as  follows,  respectively. 


[bV[d][b] 

b2\b2\  ^ 

^3i^3l 

Vi  -  ^1 

Vl  -  4'3 

bnbw  ^ 

^22^21  1 

^32^31 

Vl  - 

Vl  -</>3 

b\-ibu  ^ 

^23^2  i  1 

^33^31 

'll  -  </>: 

Vl  -  03 
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[Af]  = 


— ^  0 

!-<#>. 

0  — 

1  -  <t>2 

0  0 


where  the  transformation  matrix,  [T"],  and  its  inverse  matrix, 
[7]"'  are  given  by 


/(I  -  ^,)(1  -  <i>i) 


cos^^ 

sir\^9 

-2  cos  ^  sin 

sin^^ 

cos^B 

2  cos  9  sin  9 

«  cos  9  sin  9 

-cos  9  sin  9 

cos^B  -  sin^^  - 

[Af]  = 


M]  =  i( — L— +  — L_\ 


1  -L_.-U- 


1  -  <t>2 
0 


1  _ 1 _ 

(1-^,)=  2(1  -  -  </>:) 

_ 1 _  1 

2(1 


4Vl-</.,  1 


L  2(1  -<#.,)(  1  - 

cxp((#>i/2)  exp(^i/2)  5(exp(</>i/2)  exp(^,/2)  +  1)  2txpi^]/2) 

[M]  -  5(exp(^i/2)  exp(i>i/2)  +  1)  exp(02/2)  exp((?>:/2)  2  exp (<^2/2) 

cxp((^]/2)  cxp(<?)2/2)  2 


where  the  principa]  damage  values,  and  <^>2,  are  given  by 


cos‘^  sin^^  2  cos  ^  sin 

[7]-^=  sin^^  cos^^  -2  cos  ^  sin  ^  (50) 

-  -cos  B  sin  6  cos  9  sin  9  cos'^B  -  sin*^  - 


where  9  is  given  by 


d>]]  +  <t>2 


^4>]\  <t>22Y 


+  (/)?2  ,  (46) 


^  -  tan‘ 

2 


^  2<t>n 

^(^Ij  ~  (l>2 


respectively.  Finally,  the  damage  effect  matrix  can  be  obtained 
by  coordinate  transformation.  The  complete  set  of  transforma¬ 
tion  equations  for  stresses  in  the  principal  damage  direction 
coordinate  system  is  given  by 

[d]  =  lT]{o}.  (47) 

Similarly,  the  effective  stress  vector  in  the  principal  damage 
direction  coordinate  system  is  given  by 

[n]  =  IT]{^]  (48) 


Substituting  Eqs.  (47)  and  (48)  into  Eq.  (39),  one  obtains  the 
following  relation: 

{^]  =  [r]-‘[A/][r]{a} 

=  [A/]{a).  (52) 

The  damage  effect  matrix  [M]  is  defined  as  follow's: 

im  ^  ITV^imiT],  (53) 

Using  the  principal  damage  effect  matrix  given  by  Eq.  (41 ), 
the  damage  effect  matrix  [M]  is  given  by 


c*  2c"s‘  ^  -2r^j^  -2c^5  ^  2c5^  ^  2(c^j-cj^) 

!-(/>)  1-^2  >/(  1  -  <i>,)(  1  -*  ^2)  1  ~  4^1  1"4>2  V(  1  -  4>i)(  1  -  ^2)  1-^1  l-<^2  >/(  1  -  ii)(  1  “  4>:) 

-2r-5  5^  ^  C*  2c^5^  -2c5^  ^  2c^J  ^  -2(C^5-C5^) 

1  -  ^  1  -  <^2  V(  1  -  4>, )( 1  -  <i>: )  1-4^1  1“^2  V(1  -4),)(1  -^7)  1-^1  v'd  -4>,)(1  - 

~2r^5  2fj^  2(r^5-c5^)  -2cs^  ^  2c^s  ^  -2(c^s-cs^)  2ch^  ^  2ch^  ^  jc^  - 

1-^1  1-^2  V(  1  -  ^i)(  1  -  <^:)  1  ~  1-^2  V(I  -  <^i)(l  -  ^2)  1*"^1  1  ”  4>2  -  4>i)(1  -  <^:) 
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Conclusions 

The  fourth-order  anisotropic  damage  effect  tensor,  M,  using 
the  kinematic  measure  for  damage,  expressed  through  the  sec¬ 
ond-order  damage  tensor  is  reviewed  here  in  reference  to 
the  symmetrization  of  the  effective  stress  tensor.  This  introduces 
for  a  distinct  kinematic  measure  of  damage  which  is  complimen¬ 
tary  to  the  deformation  kinematic  measure  of  strain.  A  thermo¬ 
dynamically  consistent  evolution  equation  for  the  damage  ten¬ 
sor,  together  with  a  generalized  thermodynamic  force  conju- 
Y,  to  the  damage  tensor  is  presented  in  the  paper  by 
Voyiadjis  and  Park  ( 1995).  It  is  pointed  out  that  the  principal 
damage  tensor,  should  be  used  in  the  formulation  of  the 
anisotropic  damage  effect  tensor.  Voyiadjis  and  Venson  ( 1995) 
quantified  the  physical  values  of  the  eigenvalues,  =  1,2, 
3),  and  the  second-order  damage  tensor,  for  the  unidirec¬ 
tional  fibrous  composite  by  measuring  the  crack  density  with 
the  assumption  that  one  of  the  eigendirections  of  damage  tensor 
coincides  with  the  fiber  direction. 

The  fourth-order  anisotropic  damage  effect  tensor  in  the  prin¬ 
cipal  damage  direction  coordinate  system  is  termed  the  principal 
damage  effect  tensor.  By  coordinate  transformation,  the  fourth- 
order  anisotropic  damage  effect  tensor  is  obtained.  This  fourth- 
order  anisotropic  damage  effect  tensor  has  both  physical  sig¬ 
nificance  and  explicit  form.  It  is  therefore  not  a  mere  implicit 
mathematical  expression  to  transform  the  current  stress  tensor 
to  a  symmetric  effective  stress  tensor.  Moreover,  an  explicit 
expression  of  the  damage  effect  tensor  is  of  particular  impor¬ 
tance  in  order  to  obtain  the  constitutive  relation  in  the  damaged 
material. 

In  this  work,  tensorial  forms  are  used  for  the  derivation  of 
such  a  linear  transformation  tensor  which  is  then  converted  to 
a  matrix  form.  The  explicit  matrix  expressions  of  the  damage 
effect  tensor  are  derived  for  both  three-dimensional  and  two- 
dimensional  problems.  Especially  the  damage  effect  matrix  of 
two-dimensional  problems  can  be  expressed  by  the  damage 
variables,  4>22  and  <t>i2- 
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Damage  in  MMCs  Using  the  CMC 
Part  I:  Theoretical  Formulation 


George  Z.  Voyiadjis  and  Babur  Deliktas 
Department  of  Civil  and  Environmental  Engineering 
Louisiana  State  University 
Baton  Rouge,  LA  70803 

1  Abstract 

In  this  work  the  incorporation  of  damage  in  the  material  behavior  is  investi¬ 
gated.  Damage  is  incorporated  in  the  generalized  cells  model  (GMC)  (Paley  and 
Aboudi,  [17])  and  applied  to  metal  matrix  composites  (MMCs).  The  local  incre¬ 
mental  damage  model  of  Voyiadjis  and  Park  [26]  is  used  here  in  order  to  account  for 
damage  in  each  subcell  separately.  The  resulting  micromechanical  analysis  estab¬ 
lishes  elasto-plastic  constitutive  equations  which  govern  the  overall  behavior  of  the 
damaged  composite.  The  elasto-plastic  constitutive  model  is  first  derived  in  the 
undamaged  configuration  for  each  constituent  of  the  metal  matrix  composite.  The 
plasticity  model  used  here  is  based  on  the  existence  of  a  yield  surface  and  flow  rule. 
The  relations  are  then  transformed  for  each  constituent  to  the  damaged  configu¬ 
ration  by  applying  the  local  incremental  constituent  damage  tensors.  The  overall 
damaged  quantities  are  then  obtained  by  applying  the  local  damage  concentra¬ 
tion  factors  that  are  obtained  by  employing  the  rate  of  displacement  and  traction 
continuity  conditions  at  the  interface  between  subcells  and  between  neighboring 
repeating  cells  in  the  generalized  cells  model.  Examples  are  solved  numerically  in 
order  to  explore  the  physical  interpretation  of  the  proposed  theory  for  a  unit  cell 
composite  element. 


2  Introduction 

In  the  analysis  of  composite  materials  one  can  fallow  a  continuum  approach  or 
a  micromechanical  approach.  In  the  continuum  approach,  the  composite  mate¬ 
rial  is  treated  as  an  orthotropic  or  transversely  isotropic  medium.  The  classical 
equations  of  orthotropic  elasticity  are  used  in  the  analysis  .  This  approach  makes 
no  distinction  between  the  two  different  phases  of  the  matrix  and  the  fibers.  It 
lacks  accounting  for  the  local  effects  especially  due  to  the  interaction  between  the 
different  phases. 

For  the  last  two  decades,  researchers  have  been  using  micromechanical 
methods  in  order  to  analyze  the  multiphase  composite  medium.  Using  the 
micromechanical  approach  has  distinct  advantages  over  the  continuum  approach 


1 


ill  the  sense  that  the  local  effects  can  be  accounted  through  the  volume  average 
stress  and  the  strain  increments  in  each  of  the  phases.  This  in  turn,  is  linked 
with  the  overall  composite  behavior.  Different  micromechanical  models  employ 
different  method  of  achieving  the  local-overall  relations.  Hill  [18,  19]  employed 
the  volume  averages  of  stress  and  strain  increments  in  the  different  phases  and 
introduced  concentration  tensors  to  relate  these  volume  averages  of  the  local  fields 
to  the  overall  uniform  increments.  Dvorak  and  Bahei-El-Din  [6,  7,  8]  used  Hill’s 
technique  to  analyze  the  elasto-plastic  behavior  of  the  fiber  reinforced  composite. 
They  considered  the  matrix  to  be  elasto-plastic  while  the  fiber  is  elastic  .  In 
the  micromechanical  analysis  of  elasto-plastic  composites  Dvorak  and  Bahei-El- 
Din  [8]  identified  two  distinct  deformation  modes,  the  matrix  dominated  and  the 
fiber  dominated.  The  first  mode  is  prevalent  for  the  case  of  stiff  elastic  fibers 
while  the  second  mode  is  more  treated  as  a  general  case  of  plastic  deformation  of 
a  heterogeneous  medium. 

Recently,  Paley  and  Aboudi  [17]  developed  the  generalized  cells  model  which 
is  capable  of  predicting  the  behavior  of  metal  matrix  periodic  composites  from  the 
given  properties  of  its  constituents.  The  generalized  cell  method  is  developed  on 
the  basis  of  the  method  of  cells  model  which  was  originally  proposed  by  Aboudi  [1]. 
Its  applicability  and  reliability  in  the  sense  of  composite  properties  such  as  elastic, 
thermo-elastic,  viscoelastic  response  of  composites,  and  fatigue  failure  curves  are 
discussed  by  Aboudi  [2]. 

Kachanov  was  the  first  pioneer  who  started  the  continuum  damage  mechan¬ 
ics  .  Lemaitre  [15],  Chaboche  [3]  and  Krajdnovic  [5]  used  the  continuum  damage 
mechanics  to  analyze  different  types  of  damage  in  materials  ranging  from  brit¬ 
tle  fracture  to  ductile  failure.  Researchers  have  used  continuum  damage  model 
to  analyze  damage  in  composite  materials  by  modeling  the  composite  medium 
as  transversely  isotropic.  However,  the  proposed  continuum  approach  by  these 
researchers  made  no  distinction  regarding  the  different  phases  in  the  analysis  of 
deformation  and  damage. 

In  the  micromechanical  approach,  the  local  damage  effects  are  characterized 
separately  in  the  sense  that  the  damage  tensors  are  introduced  for  each  phase 
of  the  composite  system  (Voyiadjis  and  Kattan  [23]).  For  the  two  phase  compos¬ 
ite,  a  matrix  damage  tensor,  M"*,  is  assumed  to  reflect  all  type  of  damages  that 
the  matrix  material  undergoes  like  nucleation  and  coalescence  of  voids  and  micro¬ 
cracks,  and  a  fiber  damage  tensor  ,M^,  which  is  considered  to  reflect  all  types  of 
fiber  damages  that  fiber  materials  undergo  such  as,  the  fracture  of  fiber  (Voyiadjis 
and  Park,  [26]).  In  this  research  the  interfacial  damage  effect  is  also  expressed 
through  the  fourth  order  damage  tensorM**.  Finally  the  overall  damage  tensor  M 
is  introduced  that  accounts  of  all  those  separate  damages  of  the  matrix,  fiber,  and 
interfacial  effects.  However,  in  this  work  the  interfacial  damage  effect  is  considered 
as  a  component  of  the  subcell  itself.  This  can  be  either  matrix  or  fiber  or  any  other 
material  depending  on  which  material  occupies  the  corresponding  subcell. 
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A  thennomechanical  constitutive  theory  was  recently  proposed  by  Allen  and 
Haris  [4]  to  analyze  the  distributive  damage  in  the  elastic  composite.  In  particular, 
the  problem  of  matrix  cracking  has  been  extensively  studied  in  the  literature  (Dvo¬ 
rak  et  aJ  [10];  Dvorak  and  Laws  [9];  Laws  and  Dvorak).  Recently,  Voyiadjis  and 
Guelzim  [22]  developed  an  incremental  damage  theory  for  metal  matrix  composites 
based  on  the  modified  damage  model  outlined  by  Voyiadjis  and  Kattan  [23],  and 
Voyiadjis  and  Park  [26]. 

In  this  work,  the  incremental  damage  model  by  Voyiadjis  and  Park  [26]  is 
incorporated  into  the  generalized  cells  model  of  Paley  and  Aboudi  [17]  in  order  to 
analyze  the  damage  behavior  of  metal  matrix  composites  under  monotonic  loading 
conditions.  An  attempt  is  made  here  in  order  to  obtain  damage  parameters  for 
each  subcell  micromechanically  based  on  the  incremental  damage  model  which  is 
developed  within  the  frame  of  the  effective  stress  concept  as  presented  by  Voyiad¬ 
jis  and  Guelzim  [22].  The  subcells  (local)  damage  parameters  are  then  related  to 
the  overall  damage  variables  via  the  concentration  tensors.  These  concentration 
tensors  are  derived  in  the  damaged  configuration  in  terms  of  the  undamaged  con¬ 
centration  and  the  corresponding  incremental  damage  tensors.  A  damage  criterion 
(Voyiadjis  and  Park,  [26])  is  used  here  for  the  damage  evolution  for  each  subcell. 
The  damage  evolution  mechanism  for  each  subcell  is  considered  separately,  and  the 
extremum  principle  is  used  in  order  to  formulate  the  damage  evolution  expression. 
Finally  damaged  constitutive  relations  are  formulated  in  order.  Making  use  of  the 
micromechanical  model  (GMC)  which  allows  one  to  divide  the  repeating  volume 
element  into  many  subregions  together  with  the  incremental  damage  model,  one  is 
able  to  analyze  the  damage  at  various  locations  and  at  any  increment  of  loading. 
This  method  provides  a  computationally  efficient  approach  to  predict  the  damage 
by  using  the  GMC  for  MMCs. 

The  incremental  damage  approach  coupled  with  the  micromechanical  plas¬ 
ticity  model  is  formulated  here  in  order  to  analyze  the  damage  behavior  of  com¬ 
posites  in  the  plastic  domain  as  well  as  in  the  elastic  one.  The  fibers  in  this  work 
are  aligned  and  have  a  linear  elastic  behavior  while  the  matrix  is  considered  to  be 
elasto-plastic  material  that  obeys  the  von  Mises  yield  criterion  with  an  associated 
fiow  rule  and  a  Ziegler-Prager  kinematic  hardening  rule.  However,  the  resulting 
yield  condition  for  the  damaged  composite  system  is  a  combination  of  the  gen¬ 
eralized  Ziegler  Prager  rule  and  the  Phillips  type  rule.  The  motion  of  the  yield 
surface  is  described  by  a  kinematic  hardening  rule  that  is  a  linear  combination  of 
the  Ziegler-Prager  kinematic  hardening  rule,  and  the  Phillip’s  hardening  rule  in 
the  direction  of  loading.  In  the  numerical  simulation  of  this  work,  the  unit  cell  case 
of  the  GMC  model  is  considered  such  that  the  repeating  volume  element  consists 
of  an  elastic  fiber  region  together  with  three  elasto-plastic  matrix  domains. 


3  Theoretical  Preliminaries 

3.1  The  Generalized  Cells  Model 

The  generalized  cells  model  is  the  generalization  of  the  method  of  cells  (Aboudi,  [1]) 
by  taking  any  number  of  subcells  rather  than  four  subcells  and  considering  the  rate 
dependent  relations  of  the  subcell  for  modeling  the  multiphase  composite  materials. 
This  generalization  is  particularly  advantageous  when  dealing  with  elasto-plastic 
composites,  since  yielding  and  plastic  flow  of  the  metallic  phase  may  take  place  at 
different  locations.  The  CMC  is  able  to  provide  a  more  accurate  representation  of 
the  actual  microstructure. 

This  micromechanical  analysis,  based  on  the  theory  of  the  continuum  media 
in  which  equilibrium  is  ensured,  can  be  summarized  essentially  as  follows.  A  re¬ 
peating  volume  element  of  periodic  multiphase  composite  is  flrst  identified.  This  is 
followed  by  defining  the  macroscopic  average  stresses  and  strains  from  the  micro¬ 
scopic  ones.  Continuity  of  traction  and  displacement  rates  on  the  average  basis  are 
then  imposed  at  the  interfacesbetween  the  constituents.  The  micro  equilibrium  is 
guaranteed  by  the  assumption  that  the  velocity  vector  is  linearly  expanded  in  terms 
of  the  local  coordinates  of  the  subcell.  This  forms  the  relation  between  the  mi¬ 
croscopic  strains,  and  the  macroscopic  strains  through  the  relevant  concentration 
tensors.  In  the  final  step  the  overall  elasto-plastic  behavior  of  multiphase  inelastic 
composite  is  determined.  This  is  expressed  as  a  constitutive  relation  between  the 
average  stress,  strain,  and  plastic  strain,  in  conjunction  with  the  effective  elastic 
stif&ess  tensor  of  the  composite.  In  this  study  the  same  steps  are  followed  but 
in  addition  the  damage  mechanics  is  incorporated  by  using  the  micromechanical 
approach  in  order  to  obtain  the  damaged  response  of  each  constituent  as  well  as 
overall  instantaneous  damaged  behavior  of  the  elasto-plastic  composite. 

A  unidirectional  fibrous  composite  is  considered  here  in  the  method  of  cells. 
It  is  assumed  that  the  composite  has  a  periodic  structure  in  which  unidirectional 
fibers  are  extended  in  the  X\  direction.  This  representative  volume  element  is 
shown  in  Figure  (la).  The  representative  volume  element  (Figure  (lb))  consists 
of  NffhyN^  subcells  such  that  the  area  of  the  cross  section  of  each  subcell  is 
with  =  l...iV)s,7  =  and  each  subcell  has  its  own  local  coordinate  system 
(xi,a:2^^\a:3('^^)  with  its  origin  located  at  the  center  of  each  subcell. 

Unlike  the  method  of  cells,  in  this  work  the  instantaneous  behavior  of  the 
composite  is  considered.  The  displacement  rate  dot  denotes  time  derivative) 
is  expanded  linearly  in  terms  of  the  distance  from  the  center  of  the  subcell  (Paley 
and  Aboudi,  [17]).  This  leads  the  following  first  order  expression 

(1) 

where,  is  the  rate  tensor  of  the  displacement  components  at  the  center  of 

the  subcell,  and  are  microvariables  rates  that  characterize  the  linear 

dependence  of  the  displacement  rates  on  the  local  coordinates 


# 


# 


# 


4 


The  small  strain  rate  tensor  and  the  constitutive  law  for  the  material  that 
occupies  the  subcell  {0'y)  are  given  by  the  following  expressions  respectively 


(2) 


(3) 


The  instantaneous  stiffness  tensor,  depends  on  the  deformation  history, 

loading  path  and  applied  loading  rate.  In  this  study  for  elasto-plastic  materials* 
the  von  Mises  yield  criterion  with  an  associated  flow  rule  and  the  Ziegler-Prager 
kinematic  hardening  rule  are  used.  This  elasto-plastic  tensor  in  the  undamaged 
material  is  given  by  the  following  relation 


'pgij  ^klra 


where  Q  is  given  by 


(4) 


8a'f’>  asSf  I '  « 


(5) 


In  equation  (5),  is  the  von  Mises  yield  criterion  with  kinematic  hardening 
expressed  in  terms  of  the  backstress  tensor  The  material  parameter 
pertains  to  the  evolution  behavior  of  the  back  stress  (Voyiadjis  and  Kattan,  [23]). 
In  the  special  case  of  perfectly  elastic  materials  is  replaced  by  the  standard 
elastic  stifbiess  tensor  which  characterizes  the  behavior  of  elastic  materials 
in  the  subcells.  More  elaborate  plasticity  models  for  the  in-situ  characterization  of 
metal  matrix  composites  is  given  by  the  first  author  in  other  works(Voyiadjis  and 
Ganesh,  [24])  .  However,  in  this  work  a  simple  model  is  used. 

The  objective  of  the  work  outlined  by  Paley  and  Aboudi  [17]  is  to  solve  the 
microvariables  given  in  equation  (1).  This  equation  is  substituted  into  the  small 
strain  tensor  by  employing  the  rate  of  displacement  and  traction  continuity  con¬ 
ditions  at  the  interfaces  between  the  subcells  and  between  neighboring  repeating 
cells  in  order  to  obtain  the  relation  between  the  average  subcell  strain  rate  com¬ 
ponents  and  the  average  overall  strain  rate  components  via  strain  concentration 
tensors.  The  first  step  is  to  write  a  set  of  continuum  equations  in  terms  of  the 
microvariables.  These  interface  conditions  are  shown  in  Figures  2  and  3.  Since 
it  is  ensured  that  at  any  instant  the  component  of  displacement  rates  should  be 
continuous  at  the  interfaces,  the  following  relations  can  be  obtained  in  terms  of 
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the  micro  variable  rates  using  the  continuity  conditions  of  displacement  rates  at 
the  interfaces  between  the  subcells  and  the  neighboring  cell  and  these  relations  are 
given  by 


and 


w, 


(6) 

(7) 


All  the  field  variables  in  equations  (6)  and  (7)  are  evaluated  at  the  centerline 
for  the  subcell  and  for  the  subcell  {Pi).  As  indicated  in  Figure  2(b)  since 
the  interface  is  along  the  xz  direction  one  has  for  the  subcell  {Pi)  and  for 
the  subcell  {Pi)  and  the  interface  is  along  the  X2  direction.  This  relation  can  be 
expressed  by 


x^2^  =  x^^^  —  hp/2  or  x^2^=  x^^^  +  (8) 

and 


Xg’'^  =  x^^^  -  £y/2  or  X3^^=  x^^^  +  £tj/2  (9) 

By  employing  a  Taylor  expansion  of  field  variables  in  equation  (6)  together  with 
equation  (8)  and  omitting  second  and  higher  order  terms,  one  obtains 


= •if’’ + jAjj, 


( 


^  -(h) 

:a“^» 

dx2 


(10) 


A  similar  expression  can  obtained  by  using  equation  (7)  in  equation  (9)  for  the 
interface  conditions  along  3:3  direction  such  that 


These  equations  are  valid  in  the  equivalent  continuum  medium  in  which  the 
repeating  volume  element  can  be  defined  by  a  point  P.  This  mapping  procedure 
of  repeating  volume  elements  at  P  within  the  equivalent  homogeneous  medium 
eliminates  the  discrete  structure  of  the  composite.  Since  a  composite  is  subjected  to 
homogeneous  boundary  conditions,  the  behavior  of  all  repeating  cells  are  identical, 
and  a  uniform  field  exits  at  the  equivalent  homogeneous  medium.  The  governing 
constitutive  laws  of  this  equivalent  continuum  medium  can  be  established  by  the 
generalized  cells  model. 

Prom  equations  (10)  and  (11)  the  continuum  relations  can  be  written 

in  terms  of  the  microvariables  and  and  their  explicit  expression  can  be 


found  in  reference  [17].  The  composite  standard  average  strain  rate  is  given 
by 


,  rf0  Ny 

0=1 7=1 

(12) 

It  is  possible  to  derive  a  2{N0  -f  Ny)  +  Ny  +  1  system  of  continuum  relations 
expressed  in  terms  of  the  subcell  strain  rate  tensors  by  using  the  previous  7V^-I- 

Ny  continuum  equations  together  with  expression  (12;.  After  tedious  mathematical 
manipulations,  these  relation  can  be  given  as  follow 

€ll  -  €11 

/?,  7  =  1  •  •  • ,  iV^,  {Nff  Ny  relations) 

(13) 

1 

h2  = 

7  =  1  •  •  •  ,Ny  {Ny  relations) 

(14) 

1  A^7 

^33  —  7  /  .^7^33 

7=1 

/S  =  1  •  •  • ,  {N0  relations) 

(15) 

2?23  - 

0=1  7=1 

(one  relation) 

(16) 

2^13  —  f  ^2  ^1^13”’^^ 

7=1 

/3  =  !'•'  ,N0  {N0  relations) 

(17) 

1 

2^12  —  t^2 
^0=1 

7  =  1  •  •  •  ,Ny  {Ny  relations) 

(18) 

The  above  2(7V^  +  Ny)  -f  -1- 1  continuum  relations  are  expressed  in  matrix 

form  by  Paley  and  Aboudi  [17]  as  follows 

Aais  =  Jk 

(19) 

where  the  6-order  average  strain-rate  vector  is  of  the  form 

[f]  =  [f  11  >  ^22 )  ^33 )  ^23 »  f  13 »  ^12]  (20) 


and  the  6NfiN^  order  subcells  strain-rate  vector  is  defined  as  follows 


The  Ac  is  2{N0  +  N^j)  +  + 1  by  matrix  and  involves  the  geometrical 

properties  of  the  repeating  cells  while  J  is  a  2{N0  +  N^)  +  N^Ny  + 1  by  6  matrix  . 

One  now  needs  bN^Ny  —  2{Np  +  Ny)  —  1  continuum  relations  to  complete 
the  6  N^Ny  set  of  continuum  equations.  They  can  be  obtained  by  imposing  the 
continuity  of  the  rates  of  traction  at  the  interfaces  between  the  subcells  and  between 
neighboring  repeating  cells.  The  continuity  of  average  stress  rates  at  the  interfaces 
can  be  expressed  by  the  following  relations 

=  1,2,3  (22) 

and 


=  =  W  (23) 

One  can  express  the  average  stress  rate  in  the  subcells  in  terms  of  the  average 

strain  rate  by  using  the  constitutive  law  of  the  material  (equation  (3))  in  the 
subcells.  Using  equation  (22)  and  (23)  the  remaining  continuum  equations  which 
can  be  written  in  the  matrix  form  as  follow 


Ami,  =  0  (24) 

Am  is  bNffNy-  2{Np  +  Ny)  —  1  by  bN^Ny  matrix.  Am  involves  the  instantaneous 
properties  of  the  material  in  the  various  subcells.  The  6N0  Ny  continuum  equation 
can  be  written  in  the  following  matrix  form  by  combining  equations  (19)  and  (24) 

A^i,  =  Kk  (25) 

where  the  SN^  Ny  order  square  matrix  A^  is  given  in  the  form 


A 


m 


and  Jif  = 

0  ■ 

0 

J 

(26) 


One  can  now  solve  the  linear  system  of  equations  (25)  in  order  to  obtain  the 
following  expression 


where 


(27) 


(28) 


Ac  is  the  instantaneous  strain  concentration  tensor  that  relates  the  average  strain- 
rate  tensor  in  the  subcell  to  the  average  overall  strain-rate  tensor.  The  matrix  Ac 
can  be  partitioned  into  a  number  of  N^  Ny  by  6x6  matrices  as  shown  below 


w4c  = 


(29) 


# 
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is  the  instantaneous  strain  concentration  tensor  for  the  subcell  which  relates 
the  average  strain  rate  tensor  in  the  subcell  {P  7)  to  the  average  total  strain  rate 
tensor.  One  can  now  obtain  the  overall  effective  instantaneous  stiffness  tensor  of 
the  composite  by  using  the  strain  concentration  tensor  of  the  subcell  along  with 
its  respective  subcell  constitutive  equations  (Paley  and  Aboudi,  1992). 

3.2  Incremental  Damage  Model 

In  this  study,  the  incremental  damage  model  is  used  in  order  to  characterize  the 
damage  using  the  fourth  order  incremental  damage  tensor  where  (pj)  desig¬ 
nates  the  subceU  .  The  concept  of  effective  stress  as  generalized  by  Murakami  [16] 
is  used  here  in  order  to  introduce  the  damage  for  the  {Nff  by  N^)  constituents  of 
the  composite  system.  The  is  assumed  to  reflect  all  types  of  damages  that 
corresponding  subcells  undergo  such  as  nucleation  and  coalescence  of  voids,  and 
microcracks.  This  local  damage  response  is  linked  to  the  overall  damage  response 
of  the  composite  medium  through  the  stress  and  strain  concentration  tensors.  The 
elasto-plastic  stiffness  tensor  is  derived  for  the  damaged  composite  using  the  sub¬ 
cell  incremental  damage  tensors  in  the  generalized  cells  model,  and  the  relation 
between  the  subcell  incremental  damage  tensor  and  the  incremental  overall 
damage  tensor  m. 

Kachanov  [12]  introduced  a  simple  scalar  damage  model  for  isotropic  mate¬ 
rials  by  using  the  concept  of  the  effective  stress.  The  incremental  damage  model 
was  further  developed  subsequently  on  the  base  of  the  effective  stress  concept  for 
anisotropic  materials  by  Voyiadjis  and  Park  [26]  and  Voyiadjis  and  Guelzim  [22] 
In  its  formulation  three  configurations  are  assumed  namely  the  initial  undeformed 
and  undamaged  configuration  Co  ,  the  deformed  and  damaged  configuration  C, 
and  the  state  of  the  body  after  it  has  only  deformed  without  damage  C,  (Voyiadjis 
and  Kattan  [23]  )  as  indicated  in  Figure  4. 

By  considering  the  equality  of  forces  between  the  damaged,  C,  and  the  un¬ 
damaged  fictitious  configuration  the  following  linear  transformation  can  be 
written  between  the  Cauchy  stress  in  the  configuration  C,  and  the  effective  Cauchy 
stress  in  the  configuration  C 

dA  s=  aA  or  5=  (1  —  (30) 

where 

^  A-A  ,  , 

4’  =  -3-  (31) 

and  ^  is  a  scalar  (Kachanov  [12]).  In  the  above  equations  A  and  A  are  the 
areas  of  the  crossections  of  the  axially  loaded  bar  in  the  C,  and  C  configurations 
respectively.  The  term  ^  is  a  measure  of  damage.  The  concept  of  effective  stress 


as  generalized  by  Murakami  [16]  is  given  through  the  generalization  of  equation 
(30)  such  that 


O’  =  M  :  <T  (32) 

where  M  is  a  fourth-order  damage  effect  tensor  and  is  a  function  of  the  second 
order  symmetric  tensor  The  effective  Cauchy  stress  tensor  a  need  not  be 
symmetric  or  frame  invariant.  However,  the  symmetrized  effective  Cauchy  stress 
tensor  &  used  here  is  given  by  (Lee  et  al.  [13]) 

+  (^i!  “  (33) 

^  is  a  second  rank  tensoriaJ  generalization  of  the  scalar  function  <j>  given  by  equa¬ 
tion  (30).  The  stress  given  by  equation  (32)  is  frame  independent.  Utilizing  the 
symmetrization  procedure  outlined  by  equation  (33),  a  6  by  6  matrix  form  of  tensor 
M  is  derived  by  Voyiadjis  and  Kattan  [23].  However,  the  fourth  order  tensorial 
form  of  M  is  utilized  in  this  work. 

In  order  to  find  the  incremental  damage  tensor  m,  one  can  use  equation  (32). 
The  rate  (incremental)  expression  of  this  equation  can  be  written  as  follows 

d-  =  M  :&  +  M  :  a  (34) 

The  superposed  dot  implies  the  material  time  differentiation.  In  order  for  equation 
(34)  to  be  homogeneous  in  time  of  order  one  (i.e  stress-rate  independent)  M  should 
be  a  linear  function  of  &.  This  is  demonstrated  by  the  following  expression 

=  (35) 

Since  M  is  a  function  of  (f),  one  obtains  therefore 

Miik,  =  (36) 

0<Ppg 

Consequently  by  substituting  equations  (35)  and  (36)  into  equation  (34),  the  fol¬ 
lowing  relation  may  be  written  in  the  form 

^  (37) 

where  m  represent  the  fourth  order  incremental  damage  tensor  and  is  given  by 
Voyiadjis  and  Guelzim  [22] 


^ijkl  “  ^ijkl  ^  ^ijp^jkl^pg 

where 

^  _  dMijki  ^ 


(38) 

(39) 


The  explicit  expression  for  the  fourth  order  tensor  X  in  equation  (39)  is  given  in 
Section  4.4.  The  proposed  damage  model  was  used  successfully  for  both  monotonic 
and  cyclic  loads  Vojiadjis  and  Ganesh  [25] 


# 
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4  Theoretical  Formulation 


4.1  Basic  Assumptions 

In  this  work,  Co  denotes  the  initial  undeformed  and  undamaged  configuration  of  a 
single  laminate  while  Co^^^^  is  the  initial  undeformed  and  imdamaged  subcell  sub¬ 
configuration  of  a  single  laminate.  The  composite  material  is  assumed  to  undergo 
elasto-plastic  deformation  and  damage  due  to  the  applied  loads.  The  resulting 
overall  configuration  for  a  single  laminate  is  denoted  by  C.  Damage  is  expressed 
by  generalizing  the  concept  proposed  by  Kachanov  [12].  The  fictitious  configu¬ 
rations  is  obtained  from  by  removing  the  different  types  of  damages 
that  the  corresponding  subcell  (/?7)  has  undergone  due  to  the  applied  stresses. 
The  total  or  incremental  subcell  stress  at  configuration  C^^^  is  converted  to  the 
respective  total  or  incremental  stress  at  the  fictitious  configuration  through 
the  damage  tensor  or  respectively.  The  incremental  damage  tensor 

reflects  the  damage  related  that  subcell  only.  Following  this  local  damage 
description,  local-overall  relations  are  used  to  transfer  the  local  damage  effect  to 
the  whole  composite  system  in  configuration  C.  This  is  accomplished  through  the 
stress  and  strain  concentration  tensors  of  the  subcells. 

The  coupled  formulation  of  plastic  flow  and  damage  propagation  is  quite 
complex  due  to  the  presence  of  the  two  different  dissipative  mechanisms  that  in¬ 
fluence  each  other.  This  could  be  indicated  by  the  fact  that  the  position  of  the  slip 
planes  affect  the  orientation  of  nucleated  microcracks.  A  phenemological  model 
of  interaction  can  then  be  applied.  In  this  work  use  is  made  of  the  concept  of 
the  effective  stress  (Lemaitre)  [14].  Making  use  of  a  fictitious  undamaged  system, 
the  dissipation  energy  due  to  plastic  flow  in  this  undamaged  system  is  assumed 
to  be  equal  to  the  dissipation  energy  due  to  plastic  flow  in  the  damaged  system. 
The  damages  at  the  single  laminate  level  are  described  separately  by  the  damage 
in  the  subcells  according  to  the  material  in  the  subcells.  The  subcell  incremental 
damage  tensors,  is  better  suited  for  use  in  the  formulation  of  the  constitu¬ 

tive  equation  of  the  damaged  material  behavior  due  to  the  incremental  nature  of 
plasticity. 

In  this  work,  direct  tensor  notation  is  employed  whenever  possible.  The  ten¬ 
sors  are  denoted  by  the  bold  face  letters.  The  following  notation  wherever  possible 
for  tensor  operation  is  followed  throughout  the  paper  for  the  second-rank  tensors 
U  and  V  and  the  for  the  fourth-rank  tensors  C  and  D.  The  following  notation  is 
used  in  this  work  U  iV  ^  UijVij,  U.V  =  t/ij-V*,,  C  :  If  =  CijkiUu,  U.C  =  UijCuki 
and  C  i  Dsz  CijmnDmnki-  ^2  Slid  <^4  respectively  the  second-rank  and  fourth- 
rank  identity  tensors:  Ja  =  and  J4  =  ^{SikSji  -t-  SaSjk)  where  6ij  is  known  as 
the  Kronecker  delta.  The  transpose  and  inverse  of  tensors  are  denoted  by  the 
superscript  ”T”  and  ”-l”  respectively. 


4.2  Local-Overall  Relations  of  The  Damage  Tensors 

In  this  section  the  relations  between  incremental  damage  tensor  of  subcells 
(/?7)  and  overall  incremental  damage  tensor  m  of  the  composite  medium  are  de¬ 
rived  by  using  the  fact  that  the  average  damaged  stress  rates  <t  can  be  obtained 
as  the  average  sum  of  the  the  damaged  stress  rates  of  the  subcells  in  the 
damaged  configuration  and  is  given  by  the  following  relation 

fi=\  7=1 

In  equation  (40)  V  is  the  total  area  of  the  representative  volume  element  while 
is  the  area  of  the  individual  subcell  in  the  damaged  configuration.  Subcell 
incremental  damage  tensor  can  be  introduced  in  a  similar  form  to  equation 
(37)  such  that 


(41) 

where  encompasses  all  the  pertinent  damages  that  the  corresponding  subcell 
undergoes.  The  effective  subcell  Cauchy  stress  rate  is  related  to  the  overall 
effective  stress  rate  &  in  the  composite  through  the  stress  concentration  tensor 
as  follows 


(42) 


where  the  effective  stress  concentration  tensor  is  given  in  the  following  ex¬ 
pression  by  Paley  and  Aboudi  [17] 

^  ^(^7) .  ^(/»7) .  (43^ 

where  is  the  effective  stiffness  tensor  for  the  subcell,  is  the  undamaged 
strain  concentration  tensor  for  the  subcell  and  the  C  is  the  overall  undamaged 
effective  stiffness  tensor  for  the  composite.  One  can  solve  from  equation  (41) 
such  that 


(44) 

Making  use  of  relations  (42)  and  (44)  in  (40),  one  obtains  the  following  expression 

^=17=1 

This  equation  can  be  easily  written  in  a  similar  form  to  equation  (37)  where  m 
represents  the  overall  incremental  damage  tensor  which  reflects  all  types  of  damages 
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that  the  composite  undergoes  including  that  due  to  the  interaction  between  the 
subcells.  The  resulting  expression  is  given  by 


# 


# 


m  = 


^=1  7=1 


-1 


(46) 


This  expression  defines  the  cumulative  incremental  damage  of  the  composite  as 
a  function  of  its  subcell  components.However,  m  may  be  expressed  in  terms  of 
the  fiber  damage  m^,  the  matrix  damage  m"*,  and  the  damage  due  to  debonding 
m**(Voyiadjis  and  Park,  [26]). 


4.3  Damaged  Strain  and  Stress  Concentration  Tensors 

Concentration  tensors  do  not  remain  constant  as  the  compcKsite  undergoes 
damage.  However,  they  are  constant  in  the  undamaged  elastic  domain.  In  this 
work  undamaged  concentration  factors  are  modified  for  the  incremental  damage 
model  in  conjunction  with  the  hypothesis  of  the  equivalence  of  elastic  strain  en¬ 
ergy  [20].  The  effective  elastic  strain  and  stress  concentration  tensors  are  obtained 
by  using  the  generalized  cells  model.  The  subcell  strain  rate  tensors  can  be  related 
to  the  overall  strain  rate  tensor  in  the  following  way 


(47) 


where  fourth  order  tensor,  is  the  instantaneous  strain  concentration  tensor  for 
the  subcell  {P'{)  and  is  given  by  equation  (29).  The  undamaged  stress  concentration 
tensors  of  the  subcells  are  already  defined  in  the  previous  section  and  their 
relations  are  given  by  equations  (42)  and  (43). 

The  damaged  concentration  tensors  can  be  obtained  in  terms  of  the  un¬ 
damaged  concentration  factors  and  incremental  damage  tensors  in  connection  with 
the  elastic  energy  equivalence,  given  by 

=  dj7(/»7)  (48) 


or 

(49) 

A 

Substituting  equation  (41)  into  equation  (49),  one  obtains  the  following  relation 


(50) 
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The  above  equation  can  be  written  for  the  overall  behavior  in  similar  form  as  shown 
below 


k  =  [m]~^  :  €  (51) 

Consequently  by  combining  equations  (50)  and  (51)  with  equation  (47),  the  relation 
between  the  damaged  strain  rate  of  the  subcell  and  the  damaged  strain  rate 
€  can  be  obtained  in  the  form  shown  below 

^  ^07) .  i  (52) 

where  is  the  damaged  stress  concentration  tensor  for  subcell  {fiy)  and  its 
expression  is  given  by 


^  .  [^]-i  (53) 

Similarly  by  using  equations  (37)  and  (41)  with  (42),  the  damaged  stress  concen¬ 
tration  tensor  for  the  subcell  can  be  given  as  follows 

^(/?7)  _  jgm  .  ^  (54) 

where  is  the  damaged  stress  concentration  tensor  for  subcell  {Pj)  and  its 
expression  is  given  by 


:  m 


(55) 


4.4  Damage  Criterion 

In  order  to  study  the  evolution  of  damage  in  composite  materials,  one  first  needs 
to  investigate  the  damage  criterion.  The  anisotropic  damage  criterion  used  here  is 
expressed  in  term  of  a  tensorial  parameter  h  (Voyiadjis  and  Park,  [26]).  It  is  clear 
that  the  damage  mechanism  for  each  subcell  of  the  composite  materials  should 
be  accounted  separately  since  each  subcell  can  be  occupied  by  a  different  type  of 
material  in  addition  their  boundary  and  geometric  conditions  can  be  different  for 
each  subcell.  Therefore  one  single  damage  mechanism  cannot  be  considered  for  all 
subcells  in  the  multiphase  composite  medium.  The  anisotropic  damage  criterion 
based  on  the  Mroz  theory  [27]  is  generalized  by  Voyiadjis  and  Park  [26]  as  follows 

gi»y)  =  gi^){Y,K)  =  0  (56) 

or 

J'W  =  -1=0  (57) 
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where 


(58) 


y(^7)  j5  ^jjg  generalized  thermodynamic  force  conjugate  to  the  damage  tensor 
The  hardening  tensor  is  expressed  as  follows 


(59) 


wh6r6  t6iisors  End  ar6  defiD6d  for  orthotropic  niEtcrials  in  terms 

of  the  generalized  Lame  constant  X?'^\  and  the  material  parameters 

which  are  obtained 

by  matching  the  theory  with  the  experimental  results.  Voyiadjis  and  Park  [26] 
used  the  following  expressions  for  and 


m  _ 


0 

A2772(^)^' 

0 


(60) 


and 


0  0  X"’-” 

V,f ”  •=  0  0 

\  0  0 


(61) 


)  is  the  scalar  representing  the  total  damage  energy  and  is  given  by  the  following 
relation 

(62) 

Jo 


or 

^(St)  =  y(^7) .  (63) 

The  generalized  Lame  constants  are  defined  as  follow  (Voyiadjis  and  Park,  [26]). 

;^.(^)^^(/»7)(i_^.(S7))2  (64) 

are  the  magnitudes  of  the  effective  moduli  of  elasticity  along  the  principle 
axes  defined  along  the  direction  of  the  fiber  and  transversely  to  them.  In  order 
to  check  the  damaging  state  of  the  material,  the  following  four  steps  are  outlined 


below  by  (Stumvoll  and  Swoboda  [21]) 


A 

o 

(elastic  unloading) 

(65) 

o 

II 

(elastic  imloading) 

(66) 

9^^^  =  0, 

gy(^) 

(neutral  loading) 

(67) 

o 

II 

.  ■^(^)  Q 

5y(^7)  •  '  ^ 

(loading  from  a  damaging  state) 

(68) 

The  case  corresponding  to  loading  or  unloading  from  an  elastic  state  is  given  by 
relation  (65).  For  elastic  unloading  it  is  represented  by  relation  (66).  In  the  case 
of  neutral  loading  it  is  represented  by  relation  (67).  Finally  for  the  loading  case  it 
is  given  by  relation  (68)  from  a  damaging  state.  It  is  clear  from  the  above  outlined 
steps  that,  the  damage  criterion  =  0)  should  be  satisfied  for  the  state  of 
damage  to  occur.  As  mentioned  before  for  the  damage  evolution  of  materials  ,  dif¬ 
ferent  types  of  micro-mechanics  damage  are  considered  for  each  subcell  depending 
on  the  material  properties  within  the  subcell.  In  this  work,  for  an  elasto-plastic 
matrix,  the  subcell  is  assumed  to  undergo  ductile  damage  while  the  elastic  fiber  in 
the  subcell  undergoes  brittle  damage  and  their  total  energy  dissipation  is  different 
from  each  other.  For  the  elasto-plastic  matrix,  two  energy  dissipative  mechanisms 
of  damage  and  plasticity  are  exhibited.  Although  the  two  energy  dissipative  be¬ 
haviors  are  not  independent  from  each  other,  in  this  work  it  is  assumed  that  energy 
dissipative  due  to  plasticity  and  energy  dissipative  due  to  damage  are  not  coupled 
with  each  other.  The  power  of  dissipation  of  the  matrix  material  can  be  written 
in  the  following  form 


n(/J7)  =  -I-  (69) 

where  11^^'’'^'*  is  the  damage  dissipation  and  II^^^^  is  the  plastic  dissipation.  It  is 
assumed  that  plastic  yielding  is  independent  of  the  damage  process,  and  therefore 
the  later  term  can  be  replaced  by  its  undamaged  configuration  and  it  can  be 
expressed  as  follows 


where  the  term  is  associated  with  kinematic  hardening.  The  associ¬ 

ated  damage  dissipation  is  given  by 


(71) 


The  term  is  associated  with  isotropic  damage  hardening.  The  first 

author  has  also  introduced  [22]  a  coupled  incremental  damage  and  plasticity  theory 
for  metal  matrix  composites.  One  can  use  the  calculus  of  t^  function  of  several 
variables  in  order  to  introduce  two  Lagrange  multipliers  Aj  and  Aj  to  form 
the  function  in  the  following  form 

-  AS^^  (72) 


In  this  equation,  a)  is  the  plastic  yield  function  of  the  elasto-plastic  matrix 

and  is  the  back  stress  tensor.  The  term  is  the  damage  potential  which 
is  a  function  of  the  thermo-dynamic  force  tensor  and  the  damage  hardening 
parameter  One  can  extremize  the  function  to  solve  for  the  Lagrange 
multipliers  AJ^'^Wd  These  necessary  conditions  axe  given  in  the  following 

form  respectively 


and  gy(^)- 


(73) 


Making  use  of  equations  (70),  (71)  in  equation  (69)  and  using  the  calculus  of 
functions,  one  can  use  the  extremum  relations  in  equations  (73)  to  obtain  the 
following  expressions 


|(/87) 


djm 


=  0 


(74) 


and 


2  5y(^7) 


=  0 


(75) 


The  Lagrange  multiplier  in  equation  (74)  can  be  obtained  by  using  the  con¬ 
sistency  condition  for  the  yield  function  for  the  elasto-plastic  matrix  in  conjunction 
with  the  Ziegler-Prager  kinematic  hardening  rule.  The  corresponding  yield  func¬ 
tion  is  given  by 

p't)  ^  ha  -  :  {&  -  &)^'>)  -  (76) 


where  is  given  by 

and  is  defined  such  that 


(77) 


(78) 


In  equation  (78)  6^^’’'^  is  the  kinematic  hardening  parameter  for  the  elasto-plastic 
subcells.  The  consistency  condition  of  the  yield  function  in  equation  (76)  can  be 
written  in  the  following  form 


=  0 


(79) 


or 


*  a<i(w 


(80) 


This  condition  assures  that  in  a  plastic  loading  process  the  subsequent  stress  and 
deformation  state  remains  on  the  subsequent  yield  surface.  One  can  use  this  con¬ 
sistency  condition  together  with  the  equations  (77)  and  (78)  in  order  to  obtain 
in  the  following  form 


Al«  = 


1 


(81) 


(82) 


If  (^7) 

where  the  scalar  quantity  If  is  given  by 

jftM  =  :  («•<«  -  S<«) 

Equation  (75)  gives  the  incremental  relation  of  the  damage  variable  for  each  subcell. 
Similarly  using  the  consistency  condition  of  the  damage  potential  one  can 
obtain  the  parameter  The  corresponding  damage  consistency  relation  can 

be  given  as  follows 


(83) 


where  can  be  defined  as  a  function  of  g^^'’'^  =  k)  or  =  g^^'^^cr,  <f>,  k). 

Equation  (83)  can  be  written  as  follows 


rm 


+ 


0 


(84) 


By  substituting  equation  (63)  and  (75)  into  (84),  the  above  equation  can  be  ex- 
pressed  in  terms  of  the  parameter  Aj  where 

+  A,<«  (W  .  *<»’>  -  0  (85) 

One  can  solve  for  the  parameter  A2'  '  from  equation  (85)  such  that 


jg^  .  ^(S7) 


Sj(fiT)  .  ag(fl7)  ,  •xrid'r)  . 


(86) 


d<f> 


©y' 
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Using  relation  (86)  with  (75)  ,  the  incremental  damage  evolution  equation  for  the 
subcell  (/?7)  can  be  obtained  in  the  following  form 

^(^7)  ^  (87) 


where  is  the  fourth  order  tensor  such  that 


6o^e-() 

(^7) 


(88) 


The  thermodynamic  force  tensor  associated  with  damage  can  be  ob¬ 
tained  by  using  the  enthalpy  of  the  damaged  materials.  This  energy  equation  is 
given  by 


<f>)  =  -  ^(^^^(a)  (89) 

2 

where  is  the  specific  energy  due  to  kinematic  hardening.  is  the 

damaged  elastic  compliance  tensor  for  the  subcell.  It  can  be  expressed  in  terms  of 
the  undamaged  compliance  tensor  and  the  damage  tensor  such  that 

E-i^-r)  =  (90) 


The  thermodynamic  force  of  the  subcell  {0i)  is  given  as  the  partial  derivative 

of  enthalpy  of  the  damaged  material  equation  (89)  with  respect  to  the  second  order 
damage  tensor  “  the  following  expression 


y(/97)  _ 


(91) 


Making  use  of  equations  (89)  and  (90)  in  equation  (91),  one  can  write  the  thermo¬ 
dynamic  force  explicitly  (Voyiadjis  and  Park)  [26]  as  follows 


sc 


obcd 

i^y) 


(92) 


If  the  material  in  the  subcell  (Py)  is  elastic,  one  can  easily  see  that  the  gradual 
degradation  of  the  elastic  material  in  the  corresponding  subcell  is  caused  only 
through  damage  and  consequently  no  plastic  dissipation  occurs  in  the  material, 
that  is  n**  =  0  in  equation  (69).  A  similar  procedure  is  followed  as  outlined  before 
to  investigate  the  damage  evolution  for  elastic  materials. 
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4.5  Overall  Damaged  Stiffness  Tensor  for  the  Model 

In  this  section,  the  elasto-plastic  constitutive  model  for  the  damaged  multiphase 
composite  medium  is  obtained.  The  procedure  can  be  outlined  by  the  following 
steps.  First  one  obtains  the  subcell  (local)  damage  quantities  in  their  respective 
damaged  configuration  from  their  undamaged  relations  such  as,  stress, 
strain  concentration  tensors,  and  undamaged  effective  stiffness  of  composite.  These 
quantities  can  be  obtained  through  the  generalized  cells  model.  This  is  followed 
by  combining  the  (N^  by  N^)  subcell  constitutive  relations  by  using  equation  (40) 
in  conjunction  with  the  concentration  factors  in  the  damaged  configuration 
in  order  to  obtain  the  constitutive  relation  of  the  overall  composite  system  in  the 
damaged  configuration  C. 

One  can  start  with  by  substituting  equation  (44)  in  equation  (40).  The 
following  relation  is  then  obtained. 

^=17=1 

The  term  in  equation  (93)  is  replaced  with  the  relation  in  equation  (3), 
where  the  fourth  order  effective  tensor  in  the  effective  configuration  is  to  be 
replaced  by  the  corresponding  stiffness  tensor  depending  on  the  properties  of  the 
material  in  the  respective  subcells.  The  resulting  expression  is  written  as  follows 

■S-  =  r  E  E  (94) 

^  fi=l  7=1 

By  substituting  equation  (50)  and  (52)  into  (94),  finally  one  can  obtain  the  follow¬ 
ing  relation 

*  =  ^  ■■  i  (95) 

^  fi=l  7=1  ^  ^ 

or 

&^C:e  (96) 

where  C  represents  the  instantaneous  overall  stiffness  tensor  of  the  multiphase 
composite  medium  in  the  damaged  configuration  C,  and  is  given  by 

(97) 

^  0=1  7=1 

Prom  equation  (97),  one  concludes  that  the  overall  stiffness  tensor  in  the  damaged 
configuration  C  can  be  expressed  through  its  subcell  (local)  stiffness  tensors  and 
strain  concentration  factors  in  the  damaged  configuration 
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5  Numerical  Simulation  of  the  Model 

The  numerica]  implementation  of  the  proposed  model  is  done  for  the  special 
case  of  the  unit  cell  model.  The  applicability  of  the  incremental  damage  model  is 
assessed  herein  by  using  the  unidirectional  metal  matrix  composite  material.  The 
damaged  response  of  the  subcells  as  well  as  for  the  overall  composite  system  is 
obtained. 

The  unit  cell  model  used  here  assumes  that  the  unidirectional  array  of 
fibers(SiC)  extending  in  the  X\  direction  is  elastic  and  isotropic  while  the  ma¬ 
trix  (Ti-14Al-21Nb)  is  elasto-plastic  work-hardening  material  and  constitutes  the 
three  subcell  regioixs  around  the  fiber.  Table  1  gives  the  material  properties  of 
this  composite.  The  loading  is  assumed  applied  incrementally  along  the  fiber  di¬ 
rection  and  damage  is  checked  only  for  the  elastic  region.  Plastic  deformations 
are  ignored  in  this  example  and  is  the  topic  of  the  companion  paper  by  Voyiadjis 
and  Deliktas,  (1997)  .The  representative  unit  cell  used  here  can  be  described  us¬ 
ing  non-dimensional  quantities  and  the  subcell  volume  firactions  can  be  given  as  a 
function  of  its  non-dimensional  quantities  {h\,h2,i\Ai 

(11)  _  hill 

~  hi 

(21)  ^  ^ 

M 

These  non  dimensional  quantities  can  be  related  to  the  volume  fractions  of  the 
fiber  and  matrix  as  follows 

f  _  ^1^1  -m  _  hi£2  +  +  /t2^2  qqq\ 

^  ~  he  he  ^  ’ 

The  relations  between  hi  and  ^i  ,  h  and  e  are  known.  The  above  non-dimensional 
quantities  can  be  easily  calculated  from  the  phase  volume  fractions. 

In  this  work  for  simplicity,  the  fiber  and  the  unit  cell  are  assumed  square 
i.e  hi  =  ei  and  h  =  £.  Prom  this  assumption  one  can  find  the  non-dimensional 
quantities  in  terms  of  the  phase  volume  fractions  as  follows 

hi  =  h2  =  l-VJ  (101) 

Once  the  non-dimensional  quantities  are  determined,  the  next  step  is  to  follow 
the  procedure  outlined  in  section  (3.1)  in  order  to  obtain  the  strain  concentration 
tensor  of  the  subcells  and  the  corresponding  overall  effective  stiffness  tensor 
C  in  the  undamaged  configuration.  One  can  easily  observe  that  in  equation  (27), 
the  strain  vector  E,  is  reduced  from  by  1)  to  a  (24  by  1)  vector  form.  The 

matrix  Ac  becomes  a  (24  by  6)  matrix  and  can  be  partitioned  into  four,  (6  by  6) 


,(12)  _ 


,  C 


(22)  _ 


hi  £2 

he 

h2^2 

he 


(98) 

(99) 


matrices  where  each  one  of  them  represents  the  strain  concentration  matrix  of  the 
corresponding  subcell.  This  matrix  is  given  as  follows 


ii''”  = 


1  0  0  0  0 

■^21  -^22  -^33  0  0 

v43i  A32  A33  0  0 

0  0  0  i444  0 

0  0  0  0  .1455 

0  0  0  0  0 


0 

0 

0 

0 

0 

•<466. 


n(/J7) 


(102) 


More  elaborate  information  about  the  strain  concentration  matrix  can  be  found  in 
references  [1, 17]. 

The  damage  evolution  for  the  subcell  of  the  proposed  model  is  performed  by 
following  the  formulation  in  section  (4 .4).  The  tensorial  manipulation  is  preferred 
in  the  numerical  solution  in  order  to  get  more  elaborate  and  consistent  results.  In 
the  damage  analysis  of  materials,  the  main  objective  is  to  satisfy  the  consistency 
condition  {g  =  0)  at  any  state  of  damage.  This  phenomenon  can  be  explained  as 
follows.  Loading  of  the  material  by  an  increment  of  stress  in  the  damaged  state 
causes  the  stress  tensor  to  move  to  the  subsequent  damage  surface,  which  defines 
the  boundary  of  the  current  undamaged  region  and  can  be  expressed  by  equation 
(77).  A  this  state  g  is  only  a  function  of  the  three  variables  a,<p  and  k.  If  the 
stress  point  lies  within  the  undamaged  region,  no  damage  takes  places,  i.e.  ^  ^  0 
and  K  =  0.  On  the  other  hand  if  the  state  of  stress  at  this  point  is  increased  by 
an  increment  of  stress,  the  current  state  of  stress  will  not  be  in  equilibrium  such 
that  g{a  +  dcr^  0,  k)  >  0  which  would  mean  that  the  current  stress  point  has  left 
the  damage  surface,  which  is  impossible. 

In  order  to  bring  the  stress  point  back  on  the  damage  surface,  an  increment 
of  damage  d<f>  and  d/c  are  induced  by  equation  (63)  and  (87)  respectively.  The 
current  damage  surface  g{er  +  der,  (f>  +  d<t>,  k  +  d«)  =  0  will  then  be  satisfied. 

The  numerical  solution  investigates  the  damage  evolution  for  each  subcell 
separately  by  using  different  damage  parameters  for  different  constituents  of  the 
metal  matrix  composite.  Small  stress  increments  are  applied  along  the  fiber  di¬ 
rection.  These  damage  parameters  for  the  matrix  smd  fiber  are  given  in  Table 
2. 

In  this  work,  three  subsequent  configurations,  the  initial  undeformed/undamaged 
Co,  the  damaged  C  and  undamaged/deformed  C  are  shown  in  Figure  4.  In  the 
one  dimensional  state  of  stress,  the  relation  between  the  scalar  value  of  the  overall 
damage  and  the  subcell  damage  can  be  obtained  by  assuming  that  the  volume 
fractions  of  the  material  in  the  initial  configuration  configuration  Co  and  in  the 
damaged  configuration  C  to  be  the  same.  The  volume  fractions  for  these  three 
configurations  are  given  as  follows 

^(^7) 

*=-7- 
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where  Co  =  c  is  assumed.  One  can  express  the  total  area  of  each  configuration  as 
a  sum  of  the  areas  of  the  subcells  such  that 

A,=  ^4,  A=i:A  (104) 

0,7=1  0,7=1  ^.7=1 

Making  use  of  equation  (31)  in  (104),  one  can  obtain  the  following  relation 

(1  -4>)A='£{1-  (105) 

4,1=1 

By  dividing  both  sides  of  above  equation  by  A  and  simplifying,  the  following 
expression  is  obtained 

if  the  term  is  replaced  by  the  initial  volume  fraction  Cp,  the  above  equation 
yields  the  following  expression 

^  Cp^^'r)^{07)  (107) 

0,7=1 

The  program  output  gives  the  damage  response  of  the  material  in  each  subcell 
as  well  as  the  overall.  In  Figure  5,  different  values  for  the  parameter  v  are  used 
to  plot  the  damage  criterion,  g,  versus  the  stress  in  order  to  study  the  sensitivity 
and  robustness  of  this  parameter.  For  the  range  of  values  used  here  1.8A10~3 
to  8A10~4  the  behavior  of  the  parameter  is  quite  robust.  In  Figure  6  different 
values  for  the  parameter  v  are  used  in  order  to  plot  the  damage  variable  (p  versus 
a  for  the  subcells  (12)  and  (22).  These  subcells  are  chosen  in  order  show  how 
the  damage  can  vary  in  each  subcell  even  though  both  cells  may  have  the  same 
material  properties.  This  implies  that  the  boundary  and  geometry  conditions  are 
effective  in  analyzing  the  damage  of  the  subcells. 

In  Figures  7  and  8  the  variation  of  parameters  rj  and  (  is  studied  by  plotting 
the  damage  versus  the  stress.  It  is  observed  that  a  0.2  change  between  the  different 
values  of  r?  is  more  sensitive  to  the  damage  behavior  of  the  material  than  a  difference 
in  ^  values  of  0.05.  The  corresponding  parameters  rj  and  ^  are  evaluated  for  the 
fiber  in  subcell  (11).  It  is  observed  that  for  r)  values  between  0.1  and  0.06  and  for 
^  values  between  0.48  and  0.52,  the  material  is  quite  sensitive  to  damage,  which  is 
indicated  in  Figures  9  and  10  respectively. 

In  Figure  11  the  damage  versus  the  stress  is  plotted  for  the  different  subcells 
together  with  the  overall  damage  in  order  to  study  the  local  versus  the  overall  rela¬ 
tion.  The  model  gives  the  expected  results  such  that  the  overall  damage  behavior 
is  the  average  of  the  local  ones. 


Finally,  the  stress  strain  curves  for  the  subcells  and  for  the  overall  compos¬ 
ite  are  plotted  and  compared  with  their  undamaged  curves  in  Figure  12.  It  is  clear 
that  there  is  a  reduction  in  the  stiffnesses  of  the  material  with  an  accompanying 
non  linear  behavior  after  the  damage  is  initiated  in  the  material. 


6  Summary  and  Conclusions 

The  main  objectives  of  this  work  is  to  incorporate  damage  in  the  generalized 
cells  model  in  order  to  predict  the  elasto-plastic  damaged  behavior  of  the  metal 
matrix  composites,  and  also  to  investigate  the  applicability  and  reliability  of  both 
the  generalized  cells  model  and  the  incremental  damage  model.  The  example 
solved  in  the  previous  section  demonstrates  the  ability  to  properly  interpret  damage 
through  the  proposed  approach. 

The  proposed  study  can  be  summarized  as  follows.  The  undamaged  consti¬ 
tutive  relations  are  obtained  using  the  generalized  cells  model.  This  model  imposes 
the  continuity  of  the  displacement  and  traction  rates  at  the  interfaces  on  the  aver¬ 
age  basis.  For  elasto-plastic  materials  within  the  subcell,  a  von  Mises  yield  criterion 
with  an  associated  Ziegler  Prager  kinematic  hardening  rule  is  used  here  in  order  to 
obtain  the  undamaged  elasto-plastic  stiffness  of  the  material  in  the  corresponding 
subcell.  The  fourth  order  damage  tensor  M  and  the  incremental  damage  tensor 
m  are  introduced  for  each  subcell  separately  (rather  than  a  two  phase).  The 
undamaged  quantities  are  then  transformed  in  to  the  corresponding  damaged  ones 
by  using  the  damage  tensors. 

The  anisotropic  damage  criterion  is  used  here  and  the  damage  evolution 
of  each  subcell  is  considered  separately.  The  challenging  part  of  this  work  is  to 
evaluate  the  eighth  order  damage  tensor.  Another  important  part  in  this  work  is 
to  to  satisfy  the  consistency  condition  at  any  increment  of  loading  in  the  damage 
state.  In  this  work,  the  eighth  order  tensor  is  calculated  correctly  by  using  MAPLE. 
The  consistency  condition  is  satisfied  by  using  the  numerical  solution  procedure 
outlined  previously  Finally  the  overall  damaged  stiffiaess  tensor  for  the  model  is 
obtained  in  terms  of  its  subcell  damaged  constitutive  relations. 

The  numerical  solution  is  performed  for  the  case  of  a  unit  cell  in  the  elas¬ 
tic  domain  by  applying  a  monotonic  increment  of  load.  The  proposed  work  is 
computationally  efficient  in  predicting  the  damaged  behavior  of  the  material.  The 
applicability  and  reliability  of  the  incremental  damage  model  has  been  established. 
Using  the  generalized  cell  model  allows  one  to  predict  more  accurately  the  damage 
in  the  subregions  as  well  as  in  the  overall  composite. 

In  this  work  an  example  pertaining  to  the  elastic  analysis  of  a  single  lam¬ 
ina  is  presented.  In  Part  II  of  the  companion  paper  of  this  work  (Voyiadjis  and 
Babur,  [11]),  the  proposed  model  is  used  for  the  elasto-plastic  analysis  with  dam¬ 
age  for  both  a  single  lamina  and  laminated  plates.  The  numerical  results  are  also 


Table  1;  Material  Properties 


Matrix(Ti-14Al-21Nbj 

Modulus 

8X10^  Mpa 

Poisson’s  Ratio 

0.30 

0.22 

Initial  Volume  Fraction 

0.65 

0*35 

Table  2:  Local  Damage  Parameters 


Matrix  Damage 

Vi 

0.08 

0.06 

V2 

0.08 

0.06 

0.08 

0.06 

El 

0.55 

0.52 

6 

0.55  1 

0.52 

6 

0.55 

0.52 

0.001 

0.001 

0.0013 

0.001 

compared  with  experimental  data. 
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